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FOREWORD 


The 1955 volume of Jndustrial Mathematics, No. 6, marks a de- 
parture from the past and a step in the further development of the 
Industrial Mathematics Society. The Society has arranged with 
Wayne University Press to publish and distribute its publication. 
Also, this year’s number is in a reduced format, a format we hope 
will evoke commendation from the membership. 


The Industrial Mathematics Society is a professional organi- 
zation whose objective it is to extend the understanding and appli- 
cation of mathematics in industry. Founded early in 1949, it marked 
the first organized attempt of any group to foster a closer relation- 
ship between mathematics and industry. 


It promotes the cause of mathematics in industry through 
monthly meetings at which technical papers are presented, through 
lectures by nationally prominent scientists and engineers, through 
panel groups which concentrate in specialized fields of knowledge, 
and through publication of worthwhile papers in its annual volume. 


Membership is open to engineers, scientists, and mathemati- 
cians who are concerned with the effective use of mathematics in 
industry. Annual dues are $5.00 per year including a subscription 
to the annual volume, /ndustrial Mathematics. Application forms 
and information may be obtained by writing to the Secretary, Indus- 
trial Mathematics Society, 100 Farnsworth, Detroit 2, Michigan. 
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The Matrix For An Industrial Game 


By Robert Davies 
General Motors Corporation 


The management of a company must often decide whether to 
manufacture a product that gives considerable prestige but that can- 
not show a direct profit. The application of game theory methods is 
illustrated by an analysis of some of the factors that must be con- 
sidered in making a decision. The alternative choices, or strategies, 
available to the company can be represented by the integers from 0 
through 10, 


0 J 2 3 4 5 6 7 8 9 10 


Strategy 1 represents a nominal effort devoted to the manufacture 
of the product, and those to the right represent successively greater 
efforts, 10 representing an all-out effort that saturates the available 
market. Strategy 0 represents no effort at all to manufacture the 
product. The cost to the company of the selection of each strategy 
can be shown directly under the number of the strategy, 


0 1 2 3 4 5 6 7 8 9 10 
0 100 120 140 160 180 200 220 240 260 280 


in which some arbitrary unit is used to measure the cost. 

The cost of the initial research, development, production engi- 
neering, and tooling to manufacture the product, even cn a small 
scale, is quite high compared with the cost of making subsequent 
additions to the output. Cost increases that come with increased 
output, of course, depend on the product, but exact knowledge of the 
increases is not essential for a rough study of the game. 

Manufacturing costs can be reduced by the value of the prestige 
that the company gets for producing the product. For example, if 
the prestige is worth 200 units when the company is the only sup- 
plier, then the net costs for the various strategies are 


0 1 2 3 4 5 6 7 8 9 10 
0 -100 -80 -60 -40 -20 0 20 40 60 80 


Negative “costs”, of course, indicate a net gain. 
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The prestige is reduced if two companies supply the product, and 
so the company’s loss depends on whether there is a competitor who 
also manufactures the product. The strategies for the competing 
company can be represented by the integers from 0 through 10 ina 
vertical column: 


0 1 2 3 4 5 6 7 8 9 10 


—ESEE ——————————————— ————— — 


0 0 -100 -80 -60 -40 -20 @Q@ 20 40 60 80 
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The resulting matrix has 11* spaces, one space for each way in 
which the two companies can select one of the 11 strategies. 

The costs that were given previously are for the case in which 
no other company manufactures the product, and thus they correspond 
to the 0 strategy in the column. [If a competitor manufactures the 
product, then the cost to the original company increases because the 
decreased prestige must be shared by two companies. To make a 
cursory analysis of the game, we do not need complete details of how 
the cost increases. We can assume that the 200 units of prestige 
credit are reduced to 100 units in the case of two suppliers and that 
the 100 units are divided between the players in proportion to their 
effort. It will be convenient to call the original company Player B 
and the new competitor Player A. The complete matrix is shown in 
the first table on the next page. 

If the costs are comparable for the two players, the matrix for 
the loss to Player A can be obtained by interchanging the rows and 
columns of the previous matrix as shown in the second table. 

The effect of a third player can be determined by introducing a 
third perpendicular axis of strategies. The two square matricies 
for two players then become three cubic matricies for three players. 
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Loss for Player B 


GAME 





Player B strategy 
2 3 4 5 








“a 

0 | 0 -100 

110 50 

2/0 66 

3/0 15 

4/0 80 
0 83 

6|0 86 

7/0 88 

s|o 89 

910 90 

10 lo. 91 
“ar 

‘ol oo oO 

1/-100 50 

2 

3 

4 

5 

6 

7 

8 

9 

10 





-80 -60 -40 -20 


93 65 80 97 
70 80 93 100 
80 90 103 118 
87 97 110 124 
91 102 116 130 
95 107 120 135 
98 110 124 138 


100 113 127 142 


102 115 129 144 
103 117 132 147 


Loss for Player A 





Player B strategy 
2 3 4 5 
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The losses of one player are not equal to the gains of the other 
player; therefore the game in this form is a non-zero-sum game. 
There is no satisfactory treatment of general non-zero-sum games. 
One difficulty in analyzing such a game lies in the determination of 
the pay-off of the game, that is, in the criterion used to determine 
which player is the winner and by how much he wins. For example, 
the two players may cooperate to minimize their losses, and so they 
both play strategy 0. On the other hand, each player may want to 
maximize the losses of the other, in which case they both play strat- 
egy 10. Finally, each player may pursue an objective that is a type 
different from the other player’s objective. 

The game can be changed to a related zero-sum game by as- 
suming that each player tries to minimize the difference between 
his loss and his opponent’s loss. The matrix for this modified game 
is obtained by subtracting corresponding elements in the two previ- 
ous matricies: 


Relative Loss for Player B 





Player B strategy 





Se ae Se. | Ee tte se 

0| 0-100 -80 -60 -40 -20 0 20 40 60 80 

1/100 0 -13 -10 0 14 28 44 62 80 98 

2| 80 13 0 O 6 18 30 44 60 76 94 

3| 60 10 0 0 6 16 26 40 54 70 86 
PlayerA 4/40 0 -6 -6 ©O 8 20 32 46 62 177 
strategy «| 20 -14 -18 -16 -8 0 10 24 36 52 66 
6| 0 -28 -30 -26 -20 -10 0 12 26 40 56 

7|-20 -44 -44 -40 -32 -24 -12 0 14 28 42 

81-40 -62 -60 -54 -46 -36 -26 -14 0 14 28 

9|-60 -80 -76 -70 -62 -52 -40 -28 -14 0 14 

10 |-80 -98 -94 -86 -77 -66 -56 -42 -28 -14 0 


The maximum relative loss to Player B for each of his strategies 
can be found by taking the maximum value in each column of this 
modified game matrix. These values are 


Player B 0 i ae ae 5 6 7 8 9 10 
strategy 


Maximum loss 100 13 0 
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Player B can reduce his relative loss to 0 by playing strategies 2 
and 3. Since the game is symmetric, Player A can also reduce his 
relative loss to 0 by playing strategies 2 and 3. Each player can 
expect a pay-off of 0, which is the value of all symmetric games. 

Each entry in column 2 of the modified game matrix is as small 
as the corresponding entries in the subsequent columns. Hence, 
Player B must never play the strategies after 2 because for each of 
his opponent’s strategies, he loses less, or makes more, by playing 
strategy 2. By the same reasoning, Player A must never play a 
strategy after 2. 

The elimination of these obviously unwise strategies reduces the 
game matrix to 


Relative Loss for Player B 
Player B strategy 


0 1 2 
Player A 0 0 -100 -80 
strategy 1 100 0 -13 ° 


2 80 13 0 


The maximum values in each column are, respectively, 100, 13, 
and 0. For each of the three strategies, Player B’s loss is certainly 
as small as the corresponding maximum 100, 13, or 0 in each case; 
Player A cannot possibly make him lose more. Player B ought to 
play the pure strategy 2 in order to keep his maximum loss toa 
minimum; that is, he is assured of no loss when he plays strategy 
2. He would, of course, like to have a smaller (negative) loss which 
would indicate a gain. Likewise, the minimum values in each row 
are -100, -13, and 0, and Player A can inflict a loss at least as great 
as any one of these numbers by playing the corresponding strategy. 
Player A must play strategy 2 and make B’s minimum loss a maxi- 
mum. Thus, each player plays the pure strategy 2 and each expects 
to win 0. 

If each company wants to minimize the difference between its loss 
and its competitor’s loss, and if the cost of producing the prestige 
product is as shown in the first matrix, then each company ought to 
put a fairly nominal effort into producing the product. There is no 
need to keep secret how much research and development work is 
being done, because each company can predict what strategy the 
other company ought to play. However, if one company selects a 
nonoptimum strategy, that fact must be kept from the opponent. For 
example, if Player A selects strategy 0 in place of 2, and if Player 
B knows of this choice, then Player B can reduce his losses to -100 
by playing strategy 1. 
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The Effect of Bandpass Limitations 
On Analog Computer Solutions 


By Joachim Kaiser, Jr. and Edward J. McGlinn 
Bendix Aviation Corporation 


1, INTRODUCTION 


In recent years, the use of analog computers for the solution of 
mathematical engineering problems has become accepted practice. 
Many industries are using small analog facilities to help solve 
simple dynamic problems. One of the simplest systems of differ- 
ential equations often solved by the use of electronic analog com- 
puters can be described by expressions of the form 


(1.1) y + 2twy + w’y = Fit). 


Furthermore, many complex systems are built up by combining 
several systems of the form 1.1 into a more elaborate system hav- 
ing in general oscillatory characteristics. Thus, it is important to 
understand the limitations placed on the solutions of such systems 
by the inherent capabilities of the electronic analog computer em- 
ployed. A major limitation is imposed by the bandpass character- 
istics of the computer components, typical of which are amplifiers 
and integrators. 


2. DAMPING ERRORS 


The system 1.1 has a typical solution of the form 


’ 1 ns o- Ct ain ah 
(2.1) y F)* ov a: sinw/1-{"t, 


where ¢ is called the damping ratio, w the natural (undamped) fre- 
quency, and * the convolution operator defined by reference 1; 


ge: 
F, (t)* F2(t) J F, (t- 7) F.(r)dT. 
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FIGURE 2.1 - AN ANALOG COMPUTER “SET-UP” FOR THE SOLUTION OF 


¥+2ley+e7 y= FW 


An analog computer “set-up” for the solution of system 1.1 is shown 
in Figure 2-1. Ideally the transfer functions of the integrators are 
given in Laplace transform form by 
(2.2) -=, 

Ss 
However, due to the bandpass limitations of the integrators, the 
actual transfer functions are more closely represented by 


(2.3) - 


Similarly for the amplifier the actual transfer function is 


1 
ed “Ts +i’ 
where the time constant T is related to the bandpass of the amplifier 
by T = 1/2 7 f;, and f,, = bandpass frequency corresponding to the 
- 3 db point in c.pss. Typical values for commercial amplifiers 
range from T = 10~ to T = 10~° sec. The bandpass limitations cause 
the effective damping ¢. of the computer solution to differ from the 
ideal damping ratio ¢ in the following way: 














whe 


the 


fier 
the 
ers 
use 
the 








THE EFFECT OF BANDPASS LIMITATIONS 9 
(2.5) t. = $+ Tw (3/2 + o*). 


Equation 2.5 is derived as follows: 


The characteristic equation for the ideal system 1.1 is 
(2.6) C(s) = s+ 2tws+w’* =0, 

with roots 

(2.7) Sn=--fy tjw/i-? . 


The characteristic equation solved by the computer set-up of Figure 
2-1 is 


(2.8) €(s) = (Ts+1) (Z* + 2fwZ} +w* =0 


where Z = s (Ts+ 1). 


Equation 2.8 is derived by referring to Figure 2-1 and writing 


(2.9) Eo = Ej Ww’ /Z,) , 
but 
(2.10) E; = Ej (-2¢w/Z) - E,/(Ts+1) + f(s). 


Combining 2.9 and 2.10 and solving for E,/f(s) produces 
(2.11) Eo / f(s) = 1/ [(Ts+1)(Z?+ 2twZ)). 


The denominator of 2.11 is recongized as 2.8. 


We rewrite 2.8 as follows: 
(2.12) C(Z) = - TsZ(Z+ 2fw) 


where C (Z) is the characteristic equation of the system 1.1 with s 
replaced by Z. Assuming that the errors in the solution of 1.1 


caused by the analog computer are to be small, we expect to find 
two roots of 2.12 of the form 


(2.13) 


7) 
By 
3 


with 


(2.14) lel<<|s,|. 
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The other three roots can be shown to have large negative real 
parts, so that error terms arising from them damp out rapidly. 
(See reference 2.) We expand C (Z) about its zero Z = s,, to obtain 


(2.15) C(s_) + (Z-s,) C'(s,) + (Z-sq) C"(sq)/2 +--+ = -Ts Z(Z+22 w). 


Setting s,,' = Ss, + e and using the inequalities 


(2.16) |2Ts,, | << 
and 
(2.17) Te’| << le 
we obtain 
(2.18) e(s,+tw) + Tsp (Spt tw) = -3TSn (Sn*+2 fw). 
Finally 
e= -Ts) - 5Ts) (S,+2tw)/(s,+tw) = -Ts, (3/2+ fw/(s,+ tw). 


Using 2.7 the real part of s, is found to be 
(2.19) R(sp) = - tw+ Tw (3/2 + f°). 


Figure 2-2 is a plot of the per cent critical damping error € fcr 
obtained from 2.5 by neglecting the term t? w. In addition the figure 
shows a set of curves obtained experimentally for two types of com- 
ponents with time constants T, = 10~ sec. and T, = 1.5 x 10™ sec. 
It is noted that for loop gains w greater than 100(rad/sec)* optimum 
results are obtained for a given loop gain when the maxirmum possi- 
ble gain is obtained from the integrators. The ordinate of Figure 2 
represents the per cent critical damping, ¢., = 1. The correspond- 
ing per cent error € {ina given damping ¢ is obtained for a given 
frequency by dividing the ordinate by ¢. For example: From Figure 
2-2 we find for a loop gain w’* = 10* = 16 c.p.s. and a gain distri- 
bution Rj/a = (Integrator Gain)/(Amplifier Gain) =1, the per cent 
error in a system ideally damped at ¢= 0.1 to be 0.1/0.1 =1 per 
cent, i.e., the effective damping obtained will be é. = 0.099 for units 
with T = 1.5 x 10~ sec. 

Thus, when simulating systems with loops of high frequency and 
low damping, it is important to use amplifiers and integrators with 
good high frequency bandpass characteristics. The method outlined 
above can serve as a guide to ascertaining the expected accuracies 
of a computer solution. 
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3. ELECTROMECHANICAL SERVO LIMITATIONS 


Another area where bandpass limitations may limit the computer 
application is in the electromechanical servos most commonly used 
for multiplication and sine-cosine resolution. 

A particular example is furnished by the design criteria de- 
veloped for the Bendix Three Axis Flight Simulator. This device is 
a computer capable of simulating the motions of aircraft in flight 
and providing signals to position a Flight Table upon which motion 
sensitive instruments are mounted. 

When simulating the roll loop of an aircraft or missile by use of 
the Flight Table, the roll servo of the table becomes an integral 
part of the simulated loop, and it is important that the bandpass of 
the servo have small effects on the performance of the simulated 
loop. 

The loop to be simulated by the Flight Simulator is characterized 
by the open loop transfer function 


(3.1) Gor = Wm / 8 (8 + 20, Wm) 


where €,, is the missile damping and w,, the natural frequency. 
Actually, however, the servo transfer function enters into the com- 
putation sothat the open loop transfer function obtained by the Flight 
Simulator is given by 


(8.2) Gor =[w3,/S(St2etpn)] *{eof/(s*+2g,uys+e)}, 


where the factor in braces represents the transfer function of a 
typical electromechanical position servo with natural frequency w,; 
and damping ¢;. It is desired to obtain a criterion for the required 
bandpass of the servo in terms of the frequency ratio R =; /w,, 
and the required servo damping ¢;. 

Setting w; /w,, = R and multiplying out the denominator of 3.2 
we obtain 


D[Gor(s)] = s* (1/R*w%,) + s° (22¢/R’wh, + om/R’ win) 
(3.3) 
+s (1/wm t 2 om ee /Rwyrn) t S(gm/Wm). 


To determine the phase and gain margins we sets = jW,, re- 
group terms in 3.3 to obtain 


(3.4) D[Gor (jwy,)] = (1/R?=1-2 6, S¢/R) +5 (Spr? Sy/R- o,/R?). 


The phase margin ¢ of the transfer function Gory is plotted in 
Figures 3-1 and 3-2 for two values of ¢. = 0.1 and ¢.. = 0.2. The 
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ordinate is the frequency ratio R=wW,/w,,, the abscissa is the 
phase margin, and the servo damping ¢, is a parameter generating 
the family of curves, These results were obtained under the as- 
sumptions that in the region of interest on the Nichol’s chart (refer- 
ence 3), 


a IGo I> IGorlv!l 
b. The line a = - 90° is approximately equal to the line 
ISorl 7% 
c. The tangency of the curves Go7 (jw) and a = 90° takes place 
on the line |Go7l = 1 


The curves of Figures 3-1 and 3-2 show that the inclusion of an 
electromechanical servo in the simulation of an otherwise stable 
system can cause the simulated system to be unstable, Further- 
more, the computer user can obtain an estimate of the system 
errors caused by the bandpass and damping of the electromechanical 
servos employed in the computer. It is interesting to note that the 
lower values of servo damping ¢; < 0.5 give better results. This is 
contrary to the usual concept of optimum servo damping which is 
about ¢, = 0.7. 


4. CONCLUSIONS AND SUMMARY 


a. Amplifier and integrator high frequency bandpass limitations 
can Cause errors in the damping of second order systems simulated 
by electronic analog computers. The errors may be estimated by 
the use of charts of the type shown in Figure 2-2. 

b. Bandpass limitations of electromechanical servos used in the 
simulation of second order loops may cause a theoretically stable 
loop to become unstable. Charts similar to those of Figures 3-1 
and 3-2 enable the computer user to predict system errors and 
serve as a guide to servo bandpass and damping requirements. 
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A Procedure For Obtaining Periodic Numerical 
Solutions To Simultaneous, Non-Linear, First Order 
Differential Equations 


By Richard J. Sullivan 
General Motors Corporation 


Recently I was confronted with the task of obtaining a numerical 
integration of a set of two simultaneous differential equations. The 
equations were of the form: 


dc _ 

7 * f, (t, C, B) 
(1) 

dB _ 

dt - f, (t, Cc, B) 


The equations were non-linear in C and B, and no analytic so- 
lution was deemed possible. Normally a point-by-point numerical 
solution could theoretically be achieved with a minimum effort on an 
electronic digital computer. We need only specify initial conditions 
Co and Bo for the two dependent variables at t = O. 

The problem became interesting in that the set of equations de- 
fined only the steady-state, periodic motion of a physical system. In 
the derivation of these equations of motion, the transient damping 
and restoring forces were ignored. Consequently, a periodic solution 
to the system (period being known) could not be expected when arbi- 
trary initial conditions C, and By were employed. 

The point to be made is that initial conditions had to be found 
that caused the variables C and B to return to their initial values at 
the end of their known period. Only then could it be said a true or 
real solution had been obtained. 

This paper will attempt to describe a technique which Patrick C. 
Hayes, a fellow I.M.S. member, and I developed to obtain permissi- 
ble initial conditions for these types of equations. In practice, the 
method is a systematic adjustment of the initial conditions of the 
dependent variables until periodicity in the solutions is observed. 
Since numerical integration is used, a medium or high-speed com- 
| puter would be needed to perform the several integrations required 
at each phase of adjustment. The method was fairly successful in 
bringing our solutions of C and B as functions of time toward simul- 
taneous periodicity. 





RICHARD J. SULLIVAN 
The system (1) implies solutions of the form: 


Cc 
B 


%, (t, Co, Bo 


“ b, (t, Co Bo) . 


" 


We define Cy and Br to be the values of C and B respectively 
when t = T, at the end of the known period. Obviously, periodic so- 
lutions require that C7 = Co and By = Bo. Therefore, let us de- 
fine C = (Cy - Co) and B = (By - Bo), these quantities being the 
amounts by which the solutions fail to be periodic. Certainly we 
may expect to observe such “gaps” upon using arbitrary initial con- 
ditions. Causing these “gaps” to vanish is a statement of the prob- 
lem under consideration. 

Briefly, the technique involves a choice of arbitrary initial con- 
ditions from the allowable Co = Bo region. The functions f, and f, 
are then numerically integrated. We may then vary Co by a small 


amount A Co and integrate once again using Cg + ACo and Bo as | 


initial conditions. As a final step at this stage, we vary Bp by ABo 
and integrate using initial conditions Co and Bo + ABo. The re- 
sults of these integrations of the C and Bcuryes at the point t = T 
determine a new set of initial conditions Cg” , Bo” which, it is 
hoped, will bring the C and B curves closer to periodicity. 

As we vary one or the other of the initial conditions, equation 
set(2) indicates we may expecta different set of solutions of C and B 
over the range (t = O, t= T). Let us call the new solutions CY and 
BY where the superscript “v” means one of the initial conditions, Co 


or Bo, has-been varied to a new value. Cy} and BY will represent | 


the respective values of C’ and B” at t = T. 

With an initial condition varied, we can, in general, expect a new 
difference between the ordinates at t = QO and t = T for both the C 
and B solutions plotted in a rectangular coordinate system. This 
new “gap” for the C curve with varied initial condition can be called 
CY, in which case CY = Cl. -(C, + AC,). This ‘gap,” CY, will be 
different from the “gap” for the C curve with unvaried initial con- 
ditions. Call this change AC. 

Thus, 


A462 -T = [Ch - (Co+ AC)] 
(3) 
-([Cyp- Co] = Cp -Cy- ACo. 


Likewise, for the B curves with varied and unvaried conditions we 
can define a change in “gap” size at t = T. 
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4B2pB’-5 


Ww 
— 


B..- (B+ AB,)] 
(4) T Oo Oo 


- [By - Bo] = B, - By - ABo 


It is understood that when Bo alone is varied, ACg =O and when 
Co alone is varied, ABo= O. 

Using set (2), we now proceed to define functions over the entire 
range (t = O, t = T): 


or 
' 


C - Co = & (t, Co, Bo) - Co = &, 
B - Bo - 2 (t, Co, Bo) - Bo = >, « 


(5) 


wr 
TT 


All we have done, in effect, is to translate the origins of the C and B 
curves to the respective initial values by subtracting either the Co 
or the Bo from each function. We now recognize periodicity in set 
(2) by having ¢, (t = T) = d (t = T) = O in set (5). 

It should now be observed that AC is the change in $, at t=T 
obtained by using first one then another set of initial conditions. 
Likewise, A B is the change in $, at the end point t = T. 

Choosing proper initial values Co and Bo, which cause $, and d> 
to vanish at t = T is equivalent to finding the initial conditions which 
cause system (1) to yield periodic solutions of C and B. 

In order to investigate the effect of the initial values Co and Bo 


upon the functions d, and ¢,, we will want to approximate the par- 
tial derivatives 


To do this we consider a Taylor Series expansion of ¢, and ¢, about 
the point (Co, Bo) at t = T. Keeping in mind that at t = T, AC 


= Ad,, the change in $, and AB = Adz, the change in $2, we can 
write 


09 


(6) 
AB 


ul 


a4 ag, 
AC AB 
oF dBo Oo 


+ ++» (higher order terms) -:- 
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For a first approximation, we will consider ABo and ACo small 
enough to truncate the two series as shown. ABo and ACo are 
small but arbitrary changes in initial conditions. 

If Co is varied to Co + ACoand Bo remains unchanged (A Bo 
= 0), the system (6) readily yields approximations to 


a, oh 
OCo on 8Co 


upon observing the changes 4C in ¢, and ABin @, at t=T. Then 
varying Bo to Bo + ABo and leaving Co unvaried ( 4Co = O), we 
can approximate 


ov og, 
BBo *™ SBo 


upon integrating the equations (1) and observing the changes A C and 
ABat t= T. 
With these partial derivatives approximated at t = T, we are ina 





better position to discuss the corrections to the original Co, Bo | 


initial values necessary to cause the d, and d, functions to vanish at 
t= T. At this point, we return to the original curves [the last inte- 
gration of set (1)], and observe the corrections necessary on the C 
and B curves at t= T to cause Cy = Co and By = Bo. Naturally, 
these corrections are respectively (Co - Cy) and (Bo - Br). 

From equation set (5), we can observe that these corrections are 
(Co- Cp = -¢, (t = T) and (Bo - By) = -6, (t = T), an expected 
situation consistent with our requirement that ¢, and $, vanish at 
t= T. 

Calling (Co- Cr) = d€ and (Bo - Br) = dB, we can solve the fol- 
lowing system of two truncated series for the unknown dCo and dBo 
to obtain a new initial point (Co + dCo, Bo + dBo) and new solutions 
to (1). 





(Co - Cr) = a = B& aco + 5B, dBo + 
(7) - 2 
(Bo - Br) = dB = oe dCo + Oe dBo +**" 


In set (7), differential signs are used for the symbol “ 4” to dis- 
tineuish this set from set (6). ‘ 

With the new initial conditions C4) = Co + dCo, BY) = Bo + dBo, 
it is hoped that the next integration will come closer to periodicity 


in B and C over the range (t = O, t = T). Thus a step-by-step im- 
' 


provement is the objective. 
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This technique assumes the availability of fairly high-speed 
computers. In our case, the integrations were performed on a CPC, 
leaving the intermediate calculations and interpretations to be done 
by hand. These are quite simple, however. A high-speed stored 
program calculator would be able to perform integrations, interpret 
results, adjust initial conditions and continually improve upon the 
solution. 

The convergence of this technique has not been investigated 
rigorously. We might reasonably expect that the starting set of in- 
itial conditions need lie somewhere near the true values of the peri-- 
odic C and B at t= 0. The size of the arbitrary increments used to 
vary Co and Bo should have restrictions depending upon the differ- 
ential equations, the region of the (Co, Bo) point under consider- 
ation, etc. Generalization of this variation technique to other types 
of equations would appear to be rather straightforward. 

A geometric interpretation of what we are looking for might be 
interesting at this point. Let us consider the plane with coordinate 
axes Co, Bo. Any permissible initial point (C = Co, B= Boat t = O) 
may be located uniquely on this plane. Each point (Co, Bo) of the 
plane can be mapped into an ordinate above the plane equal to the 
value Cr of C at t = T by integration (numerical or even analytic, if 
possible). If we assume this mapping will yield a proper, analytic 
surface, we might depict it as shown in Figure 1. 


SuRFACE ~ 
Cy = (T,c,B,) PLANE eeu? 


“oe 
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PLANE OF EQUALITY 


B=e, 








Figure 2 


We are, however, only interested in those (Co, Bo) points that 
map into values of C;=Co. Hence, the intersection of the Cy - ' 
surface with the plane of equality Cy; = Co, when projected verti- (2) 
cally on to the Co- Bo planes yields a curve J,, the set of admis- 
sible initial points for a periodic C solution. j 

A similar surface of By versus (Co, Bo) and a plane of equality 
Br = Bo will yield a curve I, on the Co- Bo plane (Figure 2). [7 
is the collection of admissible (Co, Bo) points that yield a periodic 
B solution. 

Apparently, then, any intersection of J; with J, should yield an (3) 
initial condition effecting periodicity in both the C and B curves 
over the range (t = O, t = T). t (4) 

In our problem we were able to change the error (Co- C7) di 
the C curve from 0.128 at the first to only 0.020 at the fourth initial (5) 
point (cd), BS)). This required seven integrations, or about seven 
hours of CPC time. The advantages of a stored program calculator (6) 
readily become apparent to one working with this type of problem. 

The technique, in a word, is a systematizing of a relaxation method ; Si 
by using fundamental concepts of analysis. 

So long as some reasonably rapid method for obtaining numerical | (7) 
integration is available, this approach to solving steady-state peri- 
odic problems can be easily employed. 


~ 


~ 


; 


AN 


a 


An Iterative Method For Finding the Quadratic 
Factors of a Fourth-Degree Polynomial 


By Fred F. Timpner 
Chrysler Corporation 


Finding the roots of a fourth-degree polynomial that has no real 
roots is a problem that often arises in the solution of differential 
equations. If the fourth-degree polynomial can be broken into its 

' two quadratic factors, then the two pairs of complex roots may be 
| easily determined. This paper will develop an iterative method for 
___; determining the coefficients of these factors. 
} Given the fourth-degree polynomial: 


+ ar +A + ake An = O. 


It can be expressed in terms of its two quadratic factors. 


(1) X* + A,X* + A,X? + AX + Ay 
that _ (x? + BX + By x? + CX + Coa) ’ 
C: 7 ‘ 
erti- ) (2) X* + AsX* + A.X* + A\X + Ao 
= = X* + (B, + C,)X* + (B, + C,B, + C)X? 
j 
jality + (C,B, + C,B,) X + CB, . 
rt 
iodic ; Equating coefficients of like power terms give: 
id an 3) As = CG +B, , 
irves 
' (4) Ao = Bo + B,C + Co . 
r) af | 
nitial| ‘) A. = C,Bo + CoB, , 
seven } 
lator (8) Ao = CoBo . 
blem. 
ethod | Solving equations 3, 4, 5, and 6 to eliminate C, and C, gives: 
erical! (7) C.= A - Bi, 
peri- | 








(8) Gq =, 


(9) A. = By + A,B, > B,’ + Ay , 
Bo 


A;By - BoB, + AcB, 


(10) A, B , 
0 





Multiplying equations 9 and 10 by B, gives: 
(11) A.B, = B,’ + A,BoB, - B,"By + Ao , 
(12) A,B, = AsBo° ~ BoB, + AoB, 

Let W = f, (Bo, B,) = O 


(13) W = B.? + A,B.B, - BB,” - A,B, + A, , 








(14) OB, ° 2Bo ° A,B, = B,’ - A, ’ 
WwW 
as) SY - ap, - 23B, , 
_ OW ow 
(16) dw = OB, ° Bot OB, d B, 


Similarly let u =f, (B,, B,) =O 





(17) u = A,B,’ - B,’B, + A.B, - A,B, , 
(18) Ou - 2AB,- 2B,B, - A, 
OB, 0 o~™!1 , 
Ou . 2 
(19) —— By + Ay , 
- 2 ou 
(20) du = 55-9Bo + Sp @B, - 


With equations 13 - 20 it is possible to determine W, u, and their 
increments of change for any pair of values of B, and B,. It is now 
possible to use these functions in a manner similar to Newton’s 


method for functions of a single variable. 
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(From eq. 11) 


(From eq. 12) 
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QUADRATIC FACTORS OF FOURTH-DEGREE POLYNOMIAL 








| = (21) dW+W-=o0, 
| (22) du +u-=O, 
! 
| . _ Ow Ow 
} (23). W = BR aBy + SE aB, , 
_ Ou Ou 
(24) - U = OBo dBo + OB, dB, 


Solving (23) and (24) simultaneously will give values for dB», and 
dB, to correct the initial guess for B, and B,. Repeat this process 
until W and u are reduced to zero. Having determined the values of 
B, and B,, substitution into equations (7) and (8) will give the values 
of C,and C, for the other quadratic factor. 


1) Starting assumptions for B, and B, must not assume: 
1 
Bo = [ Ao }2 ’ 
B, = sAs . 
ae Ou Ou Ow Ow 
Th ] ll = = —=— = >> = ° 
ese values will give 5B. DB, AB. OB, O 


Example: 


X* + 5X°*+ 13K? + 16K + 10 = 0, 
As = 5, Ao = 13, A; = 16, Ao = 10, 








2) 
= Bo + 5B B, - BB,” - 13B, + 10, 
OW 2 
8B. ° 2B. + 5B, - B,’ - 13, 
Ow 
. at 5 - 
OB, Bo 2BoB; , 
u = 5Bo - 16 Bo - BoB, + 10 B,, 
u 
— = 10 Bo - 16 - 2 Bokh, 
Bo 
heir u 
now — = = By + 10. 
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Summary 
mn 2 ee . An 
‘ ew | ew du du 
w | ae ft 
ITERATION| B, B, | ap. | oB,| " | 9B. | 5B; | | oB | 
> + + + + + + > 
1 0 0 10 |-13 0 ee or 10 | .769 | 1.230 . 
2 .769 | 1.230 | 4.160 | - 6.825 | 1.953 | 2.224 | -10.202| 9.409 | .786| .617 
3 1.555 | 1.847| 1.259 | - 4.066 | 2.031 | 1.214 | - 6.194) 7.582 | .388| .157 , 
4 1.943 | 2.004] .182 | - 3.110 | 1.927| .262 | - 4.358 | 6.225 | .057 | -.002 sn 
5 | 2.000 | 2.002| .004| - 2.998 | 1.992} .012|- 4.008| 6 0 | -.002 
6 | 2.000/ 2.000} o | 0 
i = i 
A, = 5 = C, +B, , A, = 10 = C,.Bo , 
ses C.4¢2 , 10 = Co2, | whe 
ous 
8 ee 5 = C, inte 
ac 
pre: 
-. X* + SX* + 13X* + 16K + 10 pai: 
= (x? + B,X + Bo) (x* + GX + C) , 
X* + 5X* + 13K? + 16K + 10 ( 


= (X* + 2K + 2) (x? + 3K + 5) 
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An Equiconvergence Theorem 


By Samuel D. Conte 
Wayne University 


1. INTRODUCTION. Many problems of physical importance lead to 
a system of differential equations of the form 


U(x) - [A + qAx)] Vix) 


ul 
oO 


(1.1) 
V(x) + [A + qfx)] U(x) 


iT) 
oO 


where A is a complex aparameter and q,(x), q(x) are real, continu- 
ous functions with continuous first derivatives defined over a finite 
interval, say 0< x <1. The Dirac wave equations for a particle in 
a central field furnish an example of a problem which can be ex- 
pressed in the form (1.1). Associated with the equations (1.1) isa 
pair of boundary conditions which are usually written in the form 


U(0) cos a + V(0) sina 


iH] 
oO 


(1.2) 
U(1) cos 8 + V(1) sing = 


' 
oO 


where a and £8 are real constants. 

The values of A for which (1.1), (1.2) has a solution are called 
eigenvalues, the corresponding solutions are called eigenfunctions. 
A problem of great importance involving the eigenfunctions is this: 
Under what conditions can a function of x be expanded into a conver- 
gent series of these eigenfunctions ? 

If in (1.1) we set q,(x) = q{x) = q(x), the resulting system be- 
comes 


U-[A+q] V 


" 
oO 


(1.3) 
V+l[At+tq]) U=0. 


If (1.3) is written as a single second order equation, its solution can 
be obtained readily and is given by 
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iT] 


U(x) A cos§ +B siné , 
(1.4) 


Vix) = -A sin &é + Becosé , 
where A, B are arbitrary constants and 


g = ax+ [* a(s)ds 


Now for the system (1.3), (1.2) it is well known that an arbitrary 
function can be expanded into a uniformly convergent series of its 
eigenfunctions. We shall show in this paper that the system (1.1), 
(1.2) has solutions whose behavior resembles very closely that of 
the corresponding system (1.3), (1.2) with q(x) = (q, + 4q,)/2 whose 
exact solution is given by (1.4). In particular we shall prove that 
the expansion of a function in terms of the eigenfunctions of (1.1), 
(1.2) behaves as regards convergence in exactly the same man- 
ner as the expansion of a function in terms of the eigenfunctions of 
(1.3), (1.2). 


2. PROPERTIES OF THE EIGENVALUES AND EIGENFUNCTIONS. 
The system (1.1), (1.2) has been investigated by Hurwitz [1] and 
most of the results of this section are given there. However, they 
will be derived here in another manner to provide continuity for the 
material given in Section 3, 

Consider the system (1.1) over the finite interval a¢ x ¢b. 
By a solution of this system we shall mean a vector quantity U(x) 
= [u(x), v(x)] of two component functions. If F(x) = [ F,(x), F,(x)] 
and G(x) = [G,(x), G,{x) | are any two vector functions we define the 
Wronskian W(F,G) by the relation 


W(F,G) = F,(x) GAx) - FAx)G,(x) 


A necessary and sufficient condition that two solutions of (1.1) be 
linearly independent is that the Wronskian be non-vanishing. 

Let us write (1.1) in the form L U,=A U, where U,=[U, V] and 
L is the matrix operator 


No 
tive 


9,( 





be 
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Now let F , G be any two functions of x with continuous first deriva- 
tives defined over the finite interval a <x <b. Then we have 


“b *b “bd 
J FLGadx - J. GLF dx = -}j, ax [Fi:G.- GF] dx 


02.1) = W,(F,G) - W(F,G) . 


If F and G satisfy the same boundary conditions at a and b, then 
Wa(F , G) = WL(F , G) = 0, and (2.1) gives 


FLGdx - f° 


(2.2) f° 


GLFdx=0. 


Now let F = (x ,A) = [g(x ,A), @(x,A)], G= Ax, r’) = [G(x,Ar’), 
g,(x , ")] be solutions of (1.1) satisfying the conditions 


\ ata, ») = -cosa, g,{a,A) = sina , 
(2.3) 


@(a,A’) = - cosa, gfa,rA’) = sina. 


Then W,(F , G)=0. Also let X(x , A) = [X,(x, A), X{x, A)] be the 
solution of (1.1) satisfying 


X(b,A) =- cosf, X{b,A) = sing. 


Then W [(x, A) , X(x, A)] = w(A) is independent of x. If A is a zero 
of w(A) , then g(x, A) , X{x,A) are linearly dependent and 


X(x, A) = k gx, A) 


where k is a constant which can be neither 0 nor « because of the 
boundary conditions. Hence, 


g(b, A) = - cos B/k, g(b, A) = sin B/k . 
Similarly, if A’ is a zero of w(A) , @{x,A’) satisfies these same con- 


ditions at x = b. Hence (2.2) holds, and since L @{x, A) = A @{x, A), 
L Ax, A’) =A’ Ox, A’) , it follows that 


(A= Df” Ax, a) Ax, a) dx = = 2) L” [ale rd Gx, 2”) 


(2.5) 


+ (x, A) @(x,A’)] dx = 0. 
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This is the orthogonality property if » # X'; i.e. the integral over 
the interval (a,b) of the product of two eigenfunctions corresponding 
to two different eigenvalues vanishes. 

If A=o+it were a complex zero of w(A), then so would A’ =X 
=o- it be, and since (x, X) is the conjugate of “x, A), we have 
from (2.5) 

*b 2 2 a 
J, Wal? + lol? jax = 0, 
which is impossible unless g =0. Hence all zeros of w(A) are real. 
We return now to (2.1) and set F = g(x, A) , G= x, A’) to obtain 


(2.6) (AX - 2). ox, r) ax, 2’) dx = - gilb, A) 4b, 2’) 
+ @,(b, A) g(b, A’) . 


We already have 
Wy [¢, X] = wld) = - gb, A) cos - g,(b, A) sing , 


and if sin 8 #7 0 , (2.6) can be written 


sin pr - ay f° Ox, ) glx, »’) dx 
(2.7) 
= wir’) gb, A) - wld) 9,(b, A’) 


If A, is a double zero of w(A), then w(A,+it) = O(t?)as t——>0., 
Hence taking A =A,+it, A’ =A,-it the right hand side of (2.7) is 
0(t”) but the left hand side is O(t). This contradiction shows that the 
zeros of w(d) are all simple. If sin 8B = 0 the same result holds. 


3. AN EQUICONVERGENCE THEOREM, REGULAR CASE. In this 
section, we prove the following theorem: 


Theorem 3.1. Let {,(x), f,(x) be integrable over [0 ¢ x €1]. Then 
the series expansions of f,(x), f,(x) for the system (1.1) behave as 
regards convergence in the same manner as the expansions of [,(x), 
{,(x) for the system 








wh 
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il} 
oO 


ee [A + qx)] V 
(3.1) 


leuteaee 6, 
where q(x) = ; [qx) + qfx)] . 


Consider the system (1.1) over the interval [0 ¢< x ¢1] and let 
Ax, A) = [P, |, K(x, A) = [X,, X,] be solutions of (1.1) satisfying 
respectively the following conditions 


Nag = - sina, 9,0) = cosa , 
(3.2) 


lxsy = - sinf , X{1) cosp . 


Now (x, A) satisfies the first of the following boundary conditions, 


u(0) cos a + v(0) sina = 0 , 
(3.3) 


u(1) cos 8 + v(1) sinfB = 0, 
and it will also satisfy the second of these if we require that 
(3.4) y,(1, A) cos 8 + gf{1,A) sing =0. 


Now (3.4) is simply the condition that the Wronskian W(¢, X) vanish, 
i.e., 


Wy, X) = 91, A) XA1, A) - gA1, A) X(1, A) , 


y,(1, 4) cos 8 + @(1, A) sinB = wiA) = 0. 


The following asymptoptic solutions for large complex A are de- 
rived (reference 2): 


(ax A) = sin (E - a) + Oet*/jal) , 
(3.5) : 


PAX, A) 


cos (— - a) + Oetx/|al) , 











( As A) sin ( - 0, - 8) + Ofe'*//Al) , 


(3.6) « 
xs A) = cos(— - 0, - 6) + Het*/|al) , 


where 
+ 4 (ads) + ads)] ds , 


1 ;/x 
= & =Ax+5 f° (a, +4) ds. 
Let us define a function ¢ (x, A) = [¢,, ¢,], 


$, (x, A) = ee v [ oy,a) fy) + 9,ly,A) £fy)] dy 


QAx,A) fF 
+ oo { [X,f,+X,f,] dy , 
(3.7) § 


' G , (x,A) = Xa {* [ 9, f, + Q, f, | dy 


Pal x,A) I 


way x [f+ %,f,) dy, 


where w(A) = W(g, X). One may verify by direct substitution that 


$, , &, satisfy the inhomogeneous equations 


$ - [A + q,] $d, = f(x) ’ 
$, a [A . q,] $, 7 f(x) ’ 


which suggests that an expansion for f,(x), f(x) may be obtained 


from (3.7). For real A we have 


(3.8) Wig, X) = wid) = sin a[cos 9 cos 0, - sin 8 sin@,|- cos a 


[-sin 8 cos @, - cos B sin @,]| + a>) : 


= sin (8, + B - a) + >) , 
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Co 
(3. 


(3. 
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Now Hurwitz (reference 1) has shown that for real A the zeros of 
w(A) are A, =O+nT + a) where 9 = ; i. (q, + q,)ds + B - a. 


For the case sin a #0, sin 8 40 we have 


1 1 L 
(3.9) wa) sin(@, + B - a) + O1/d) ’ 


2 1 
~ sin(@, + B - a) [1 + O1/r)] ’ 


. : [1 + O(1/a)] . 


sinf®, + B - a) 


Consider the first component of equation (3.7), and use (3.5) and 
| (3.6). This gives for complex A 


X(x,A) v(y,A) fy) = {sinlé, - 0, - A) sine, - a) 
+ Oe *¥-*) 1a |)} fy) 
= {} (cost, - fy - & - B+ a) 
2 x y 1 
-cos(é, + by o - - B-a) | + Oet(+¥-*) /|, |} 


f(y) , 


at 
= {F,(x, y) + He''*¥"*) /|n|)} fy) 
where the order term includes all cross product terms. Making use 
of (3.9) we have 
ed (3.10) Xx, ») gly, A) _ Fx y)  , Mett¥-*) /Ja |) 
wir sin(@,+8-a) sin(Q,+8- a) , 
. (3.10) . Fixy)  , Mell*7-*) jr) ; 


sin@,+8-a) ~ et) 


. K(x, y) -x 
* sinld,ee-a) * Mer-*//lal) 
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The same procedure can be used to obtain similar order properties 
for X,(x,A)9,(y,A), 9,(x,A) X,(y,A), 9,(x%,A) X{y,A). We now investi- 
gate the first of the terms in (3.7): 


X,(x,A) yx 


(3.11) wa) é [ ofy,At ly) + vAy,AtAy)] dy 


. s* F (x,y)fi(y) + FAx,y)f, (y) 4 
0 sin(0,+8-a@ 


ots 4* e et(y-x) [|f,(y) | + | ffy)|] dy} . 


For large real A we have A, =@®+nza. We will now show that the 
integral of the second term vanishes over the squares in the figure 
which are spaced between the poles of w(A). 








t 
Tisonnibl -- ((n+2) 7 +0, (n+4)s +9 | 
. [ 2 2 
A=or+it 
0 
r 
SE © yo] 











D A 


If 0O< 6<x and { [\f, | + | £,|] ds exists and if we designate 


the second term in (3.11) by Kx, A), we can write 


-$ * ‘ 
0 {<— ' (if, |+ | f,1] dy} 


" 


(3.12) I(x,A) 


ot Y(t lel til ay}, 


: ry: b+ of J%, fit, 1+ 1f,1) dy} . 








wh 


Th 
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We now integrate over the square contour J" 


1 x 


r 0 in fi, (if + |f,|] dy } da 

-.% 
= 016 (itl + ltl] dy}, 
where we have used the fact that over the contour/]”, f aT = 0(1). 
This last term can be made arbitrarily small by choice of 6 and 
having fixed § the first term in (3.12) tends to 0 as n —»oo since 
e 6t . e ét 


1 
Jr 01 r ial 


2ri “[ ay, 


~ ott ot 2s 
AB or CD (n+5) 1+0 


-6 [(n+5) 7+0] 


. of 5 e Sf A ’ 
. BC or DA (n+5) 7+0 
-(n+3) 7. §¢ 
= 0 rh. S e dt} 
- -(n+5) 1-0 


-6[(n+5) 7+0] 


-(n+=) 1-0 
+ os ls e do , 
jn -(n+5) 1+0 
1 ~b[(n+5) +0] 
= 0 ing }+O0f{e rs 


A similar argument applied to the second term in (3.7) yields 


Ay s) [ XAy,A)fi(y) + Xdy,A)f(y) ] dy 


= 5! G,(x, y)fi(y) + GAx,y)fAy) 4 
sin(@,+8- a) 
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1 ‘ " 
+ Ofnr J ety) []L| + [LO] ay}, 


and again the integral of the second term on the right vanishes over 
the squares as n——->-0. Thus we have left 


— x F (x,y)f,(y) + F,(x,y)f, (y) 
Sp (xara = cfr {sf walt 


rT 27i 0 sin(O +B-@ 


(3.13) - Gxy)f, + Gf, , 
x sin(0,+B-a) 


y } da + A) 


The second component of (x, A) yields a similar term. The func- 
tions F,(x, y), F,(x, y), GAx, y), G(x, y) are analogous to those ob- 
tained in (3.10). These can be written 


1 
Fi(x, y) = 5 {cos y, - cos y,} , 
rr ; 
F,(x, y) = 2 { sin "% *+ Sin Ye} ’ 
1 
G(x, y) - 2 {cos 3 cos Yo} ’ 
1 , . 
Gx, y) = 2 {sin Ys + Sin Yo} , 
where 
n=§& -&§ -Q-B+a, 


$ 4 
% = & +8 -@-B-a, 
é 


% * oy eh * 8-8 


The integral part of (3.13) is precisely what one obtains when 
the above procedure is applied to the solution of the special system 
(3.1). Indeed the exact solution of the system (3.1) is given by 

U 


Vv 


Acos&—+Bsiné , 
-Asin § + Bcecosé , 


where — =Ax + ; § | q,(s) + qAs)| ds. Let (x, A) be the function 
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corresponding to (3.7) but applying to the system (3.1). Then we 
have shown that 


Ti rp $(%) w = icity Go(x, A) a + O(1) 
n n 
where o(1) is a function of (x, y, A) which vanishes as 4 >on 
the squares J[°,. Since the squares include the poles of the func- 
tions ¢, ¢ , the residues on the left side must equal the residues on 
the right side. Hence the expansions for f,(x), f,(x) in a series of 
the eigenfunctions of the system (1.1), has the same convergence 
properties as the expansions of f,(x), f,(x) into a series of the eigen- 
functions of the system (3.1). 
The expansions for f,(x), f{,(x) are easily shown to be 


2 
f(x) = ,_2 en Ualx) , 
(3.14) 
f(x) = a Cc, V,(x) , 
where 


cn = 5" [f, (yUAy) + £Ay)Valy)] dy, 


and U,, V, are solutions of the system (1.1), (1.2) corresponding to 
the eigenvalue A,. The equiconvergence theorem proved above then 
shows that the series in (3.14) converge uniformly provided that 
f,(x), {,(x) have continuous second derivatives and satisfy the bound- 
ary conditions 


{,(0) cos a + £,(0) sina = 0 , 


f,(1) cos B + £,(1) sin 6 SB 


l 


The uniform convergence of the series (3.14) is shown directly 
in reference 1 in another manner. 


SAMUEL D,. CONTE 
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The Mechanics of Hydraulic Brake Systems 


By W. A. Van Wicklin 
Ford Motor Company 


The purpose of this paper is to suggest to mechanical engineers 
a method of dynamic analysis for machines. This method permits 
the engineer to obtain a graphical picture of a system of dynamic 
equations for the machine elements. Several types of automotive 
braking hydro-mechanical control systems and two types of power 
brake controls will be discussed. 

A typical passenger car hydraulic brake assembly is shown 
schematically in Figure 1. A pendent-type pedal is suspended under 
the instrument panel. The pedal return spring is attached directly 
to the pedal. Between the brake pedal and the master cylinder 
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assembly is a short push rod which abuts against the piston within retu 
the master cylinder. Steel tubing connects the master cylinder with mult 
each of the four wheel cylinders located inside the brake drums on on | 
each wheel. Pistons in the wheel cylinders expand the brake shoes area 
outward against the brake drum. Friction between the drum and four 
shoes absorbs the kinetic energy of the moving vehicle and this the | 
energy is dissipated as heat. | fore 

A block diagram method of notation has been developed which is pres 
able to show the force, position, velocity, and acceleration for each beer 
member of a system. While this method is more elaborate than | plex 
necessary for the mathematical analyst, it is becoming accepted by grea 
mechanical engineers because of the ease of observing any desired amp 
variable or the effect of changing any constant. by a 

A block diagram of this type for the hydraulic brake system is » the 
shown in Figure 2. To operate the brakes, the driver applies a not | 
force, F;, at the brake pedal. To move the pedal, this input force I 
must be greater than the preload of the return spring. Therefore, ' rela 
in the diagram, the first box shows the comparison of the input and load 


STANDARD HYDRAULIC MANUAL BRAKE SYSTEM 
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return spring forces. The difference between these forces is then 
multiplied by the ratio of pedal to push radii to give the force, Fp, 
on the master cylinder piston. Dividing this force by the piston 
area, Ayyc, the hydraulic pressure, P)yc, is obtained. Since the 
four wheel cylinders act in parallel, it is convenient to represent 
the load by one wheel assembly which is equivalent to the four. The 
force, F,, on this wheel cylinder is the product of the hydraulic 
pressure and the piston area, A,. While the force, F,, could have 
been represented by multiplying F 4 by a constant, the addition com- 
plexity resulting from the previous development is justified by 
greater identification that can be made for each e‘ement. As an ex- 
ample, the hydraulic pressure would not be showi:, if multiplication 
by a constant had been used. Therefore, if our interest had been in 
the variations on pressure, P,,-(t), as F(t) was changed, it could 
not be seen in the diagram. 

Experimental static measurements have been made to obtain the 
relationship between output force and position. It was found that the 
load could be represented by a mass and a non-linear spring. The 
non-linearities of the spring take into account the compressibility of 
the brake fluid and the elasticity of the steel tubing as well as the 
return spring on the brake shoes and the deflection of the drum and 
shoe assembly. The mass includes all the moving parts of the sys- 
tem: the pistons, shoes, fluid, and pedal. 

The output position, X, is found by first subtracting the spring 
force, F,,, from the output force, F,., dividing by the mass and then 
integrating twice with respect to time. The brake pedal position is 
the product of the output position with the ratio of the wheel cylinder 
to master cylinder areas and the ratio of pedal to push rod radii. 
The pedal return spring force is the product of the pedal position, 
Xa, and the return spring rate K.. This spring force completes the 
diagram. An equation can be easily written directly from the dia- 
gram describing this system. 


(1) F; = K,(a/b)(A</Ay4OXo + (b/a)(Ayac/As)(MoXo + Ks Xo) 


Due to the non-linear term in this equation, it is convenient to use 
machine computing techniques if complete solutions are required. 
The natural frequency of the system can be found for any load. 

The most important part of the system shown in Figure 1 is the 
driver. The driver is a sensing device, controller and actuator that 
produces the input force to the braking system. The driver requires 
some positive measure of the brake rate other than by visual feed- 
back and body forces caused by the deceleration. Usually this is in 
the form of a pedal load which is a function of the pedal movement. 


42 W. A. VAN WICKLIN 


This pedal load function is inherent if the brakes are manually 
operated. 

In the power-controlled systems, there are several means of 
providing the pedal load. One method is to furnish an artificial load 
by means of springs so that the pedai load depends solely upon the 
pedal position. The only effort the driver applies is a comparatively 
light force to overcome friction plus the spring load. Such a control 
is called a Full Power System. 

Frequently in hydraulic controls, reactive-type valves are used, 
A reactive-type valve is a directional-control valve in which the 
controlled pressure to the actuator also acts on a reaction area of 
the valve which produces a force opposed to the input force and 
tends to close the valve. When this reaction force is not added to 
the power force as part of the output, the system remains a Full 
Power System. 

The other arrangement is to feed back a portion of the output 
load to the driver. This type of a powered-control is called a 
Booster System, since the driver must contribute a definite part of 
the output work, the remainder being furnished by the hydraulic 
power source. Sometimes these systems are called Hydraulic- 
Assisted. 

Booster systems are of two types: force and position feed-back. 
The pressure booster is typical of the force type. In this system the 
booster is placed in the brake line between a standard master cylin- 
der and the wheel cylinders. The booster acts as a constant-gain 
pressure amplifier; that is, the downstream brake line pressure is 
a linear function of the master cylinder pressure. 

Figure 3 shows the general physical design of a typical hydraulic 
power brake booster. Master cylinder pressure acts upon area Av, , 
on the control valve, and on area Ajyy, on the output piston. Opposing 
the hydrostatic force on the valve spool is the mass of the valve, the 
return spring, and the hydraulic pressure P;,; controlled by this 
valve, which acts upon area App. The valve position, which is de- 
termined by the forces acting on it, regulates the flow of oil in and 
out of the booster cylinder. The control valve is a reactive-type 
3-way valve. Oil is supplied to the valve at pressure P,. Pressure 
Py also acts on the area App of the booster piston and produces a 
power force which is added to the manual force on the output piston, 
The sum of these two forces divided by the area of the output piston 
gives the brake line pressure. The remainder of this system is 
identical to the manual brake system. 

This design is both a pressure amplifier and a volume amplifier; 
that is, the output volume of the booster is greater than the output 
volume of the master cylinder. This additional volume is obtained 
by the volume displaced by the booster piston rod area. 
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HYDRAULIC POWER BRAKE BOOSTER 


__ MASTER CYLINDER 
~~ CONNECTION 


BRAKE LINE 
bao CONNECTION 


sai 





Fig. 3. 


Boosters of this type are currently used on passenger cars and 
trucks which utilize engine manifold vacuum or compressed air as 
the power source rather than oil hydraulics. 

Figure 4 is a diagram of this power booster in a brake system. 
The booster part is enclosed by a broken line. Since there is no 
position feed-back of the output position to the control valve, this 
booster is not a servomechanism. Output velocity, however, is fed 
back and this provides most of the damping. Automatic feed-back to 
the control valve comes from both the valve position and pressure 
Py, which are compared with the master cylinder pressure. 

The steady-state input force and output position relationship is: 


(2) F. = AX, +BX,+C, 
where a - , tb/a)(Amc/App)Ks 





[AV,/AV2 + Am/ABp |] 
B = (a/b)(Am/Amc)(As/Ac)Ki ’ 


c - 2/a)(AMc/AV2) LKe 
~ [AV,/AV. + Aw/ABP] 
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HYDRAULIC POWER BOOSTER ADDED TO 
STANDARD MANUAL BRAKE SYSTEM 


[POWER BOOSTER = oo 


























Fig. 4. 


and L = underlap of return port on control valve 


when 
F, > K, (a/b)\(A2 Ke) /(AyycAViKs) + (b/a)(Ayyc/AV:)KeL. 


1 
If F; is less than the above expression, the system is the same as a 
manually-operated brake. Consequently, there exists a cut-in point 
for a booster system of this type. 

Dynamic equations of this system are easily written. Due to the 
pressure, flow rate, vaive position relationship, the result is a very 
non-linear pair of simultaneous, ordinary differential equations. 
Again computer solutions are best sought. 

The other type of booster system has position feed-back. An ex- 
ample is one which incorporates a follow valve or servomotor. A 
hydraulic power brake of this design is shown in Figure 5. Power is 
added to the manually applied input force before the master cylinder 
piston. This booster acts as a constant-gain force amplifier. 
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The brake pedal is attached to a reactive-type control valve | 


spool by the push rod. The input applied at the brake pedal is 
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HYDRAULIC POWER SERVO BRAKE 


f 
RESERVOIR 


BRAKE LINE. 





Fig. 5. 


opposed by the pedal return spring, the mass of the pedal, push rod 
and valve spool, and by the hydraulic reaction force on the valve 
spool, The control valve is a 3-way valve with connections to the oil 
pressure source, the oil reservoir, and the booster cylinder of area 
Ap. The booster cylinder pressure, P;,, which is controlled by the 
valve, also acts upon the valve spool reaction area Ap. Since the 
booster cylinder area and the valve spool reaction area are con- 
stants, the booster force gain is constant and equals Ap divided by 
Ap. The master cylinder pressure is then the product of the booster 
cylinder pressure and the ratio of booster piston to master cylinder 
piston areas. This pressure is transmitted through the brake lines 
to the wheel cylinder pistons. The output load is the same as in the 
other systems. Wheel cylinder piston position is fed back to the 
master cylinder piston by the column of fluid in the brake line. 
Therefcre, the master cylinder piston position corresponds to the 
output position of the system. The master cylinder piston is also 
the control valve body. Consequently, as the wheel piston displace- 
ment increases, the control valve is closed and causes the piston to 
Stop. If the input force is decreased, the valve opens to connect the 
booster piston area with the oil reservoir and reduces the pressure 
at the booster cylinder. In turn, the master cylinder pressure is 
reduced and the wheel cylinder piston displacement decreases. 
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Position feed-back of the wheel cylinder piston to the valve body 
closes off the return port to stop the motion. A new steady-state 
output position is then established that corresponds to the new input 


force. The control system error, that is, the difference between the | 


input command and the output position, is the relative displacement 





between the valve spool and the valve body. An advantage of a hy- | 


draulic control of this type is that a very stiff system is obtained, | 


Stiffness is defined as the force required to produce an increment of 
error. With an ideal, zero-lapped valve and an incompressible po- 
sition feed-back loop, the stiffness is infinite. Usually, however, to 
reduce the quiescent leakage rate, the valve is made with a slight 
overlap which produces a dead band. 

A diagram of this system is shown in Figure 6. This is a servo- 
mechanism since the automatic control is operated by position feed- 
back. If the reaction area Ap, was zero, this would be a Full Power 
System. But when Ap is greater than zero, this isa Booster System, 
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The output power of this system is H = Qp* P,, . 

But Qp = Qp +Qy, and Qp = Ze* de 

Therefore, H = Xp Ap Py + Qy Py 


Xn Ap Py is the manual portion of the output and Q,, P,;, is the 
power added by the power source. 


A set of dynamic equations which describes this system is easily 


| written: 
(3) F; = (a/b)(Ayyc/A,) KiXo + (b/a)(Pyy Ap + Mp Xp) 
(4) MoXo + KsXo = Pyy(Ap/Aysc)A, 
F.. Xp > Xp => Qy PO, 


(5a) Ap(A;/Amc)Xo - ApXp = C (Xp - AyycXo/A,) VP.-Py;. 


| or when Xp < Xp - > Qu He FF 


| (5b) Ap(A,/Ay4c)Xo - ApXp = Cy (Kp - AyycXo/As) VPy « 


With the above three equations having the input force F; , the valve 
position Xp, the control pressure P};, and the output position X,, as 
the four time dependent variables, the equations of motion have been 
obtained. The steady-state equation for input force and output po- 
sition has the same form as equation 2. 
| These power brake systems, as examples, show how hydro- 
| mechanical systems can be described mathematically. By the use of 
a block diagram, an entire system is expressed using the equations 
of dynamics, kinematics, elasticity, and fluid mechanics. One of the 
principal advantages of this form of representation is that it is an 
| aid to a better understanding of the whole system. These diagrams 
should not be confused with other similar block diagrams. The 
presence of many loops in these diagrams does not imply that the 
|}machine is a multi-loop automatic-controller. Instead, each loop 
shows the path that some force, pressure, position, velocity or ac- 
celeration follows through the machine. These loops show how the 
machine elements interact upon each other. 
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Election Forecasting on UDEC * 


By Saul Rosen 
Wayne University 


Several well-publicized attempts have been made to use elec- 
tronic computers to predict final nationwide election results on the 
basis of early returns. At Wayne University for the election in No- 
vember 1954 the Computation Laboratory made what we believe to 
be the first attempt at electronic computer forecasting of the results 
of a statewide election, using a relatively large sample and attaining 
extremely accurate results in a very close election, 

The project was carried out in cooperation with reporters Lou 
Arkles and Al Creecy and other members of the staff of the Detroit 
Times. They arranged to collect election returns, when available, 
from precincts scattered all over the state of Michigan, and set up a 
direct phone link with the Computation Laboratory to feed in the data 
as soon as it was available. 

The state was divided into four sections: Detroit, the tri-county 
metropolitan area (Wayne, Macomb, Oakland) outside Detroit, south- 
ern Michigan excluding the above three counties, and northern 
Michigan. The division between north and south was the traditional 
Muskegon-Bay City line. At first the Northern Peninsula was to be 
a fifth section, but this was not feasible since early returns from the 
far north would be too sparse, and the area was therefore included 
' innorthern Michigan. 

In preparation for the attempted forecasting the results of the 
previous two comparable elections were analyzed. For each pre- 
cinct the Democratic percentage of the vote was computed, and from 
this figure was subtracted the Democratic percentage for the sec- 
tion as a whole. This was done for both previous elections and the 
average of the two numbers thus obtained was taken as the differ- 
ential figure for that precinct with respect to its section. The dif- 
ferential is a measure based on past performance of the bias in 
| fevor of the Democrats that would exist if the results in that par- 


} ticular precinct were used to predict the results for the whole 
section. 





*Unitized Digital Electronic Computer. A large scale general 
purpose electronic computer built by the Burroughs Corporation for 
the Wayne University Computation Laboratory. The Burroughs Cor- 
poration provided financial assistance toward the election forecast- 
ing program. 
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On election night, returns from precincts were accepted as they 
came in. In each section a weighted average of the differential was 
computed, the weights being the total votes counted in each precinct. 
This weighted average of the differentials was a measure of the bias 
of the sample up to that point. It was subtracted from the Demo- 
cratic percentage in the section at that point to give a current pre- 
diction of the Democratic percentage. 

In order to get a statewide prediction it was necessary to predict 
the total vote in each section. This was done by multiplying the total 
vote in a section in 1950 by an expansion factor equal to the total 
vote in the precincts thus far reported divided by the total vote in 
the same precincts in 1950. Using the predicted total votes and the 
predicted Democratic percentage figures a statewide prediction 
could be made. 

In more concise form the following are the calculations that were 
performed: 


Previously stored in the computer were Vg, the total vote in the 
section in 1950 for governor and for each precinct: 
Vio = total precinct vote for governor in 1950, 


d; = differential in gubernatorial percentage 


as explained above. 


1 


As each precinct reported there were fed into the machine: 


v;, = total precinct vote for governor, 
% = total precinct vote for senator, 
w. = total precinct vote for Williams, 


= total precinct vote for McNamara. 


Williams and McNamara were the Democratic candidates for 
governor and senator respectively. Votes for Republican candidates 
are, of course, implicit in the data, and Republican percentages were 
finally found by subtraction from 100. 

As each precinct reported, the computer determined which sec- 
tion it belonged to on the basis of its code number, and for that 
section it accumulated the values of: 


PW ~YV: ~ Ww. > . ‘ ” . 
LVin» UVig UW, UM; TVigg ZV;,4; , 
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ELECTION FORECASTING ON UDEC 


th ZW; “7 is the new predicted Democratic 
os ZV, percentage for governor, 
= mj F is the new predicted Democratic 
Z Vis percentage for senator, 
Sy. 
Vp: = Vso se is the predicted total vote for governor, 
4 “—e (2 Viz is the predicted total vote for senator, 
- pe xv P! all forthe section under consideration. 


Using the Vp as weights the percentages for the four sections 
were combined and a statewide prediction obtained. 

The program was set up so that a new sectional and statewide 
prediction could be obtained with each new precinct, if desired. 
Thus the convergence of the predictions could be monitored. Al- 
though no formal statistical estimates of validity were attempted, 
confidence in the results increased as the estimates became more 
and more stable and converged toward the final values. 

The selection of the precincts for the sample used was not nearly 
as much under our control as we might have wished. It was neces- 
sary to use precincts for which there was a reasonable chance that 
early returns would be available. Even so there was no guarantee 
that the returns actually would be available. However, it was pos- 
sible to eliminate precincts for which past behavior was erratic, and 
also precincts which had undergone major population changes or 
even geographical changes. 

The past history for 1,000 precincts was stored in the computer. 
This is about one-fourth of the total in the state. Actually, about 400 
precincts were used when the final prediction was made at about 
1:00 a. m. Wednesday morning. Based on about a 10 per cent sam- 
ple, mostly from urban areas, UDEC predicted that Williams would 
have 55.78 per cent of the total vote in the state. He actually had 
55.76 per cent. UDEC predicted that McNamara would have 51.08 per 
cent of the total vote. He actually had 50.92 per cent. The UDEC 
predictions were made at a time when McNamara was far behind on 
the basis of current tabulations, and publication of the prediction in 
their first edition represented a major scoop for the Detroit Times. 

The method used was quite elementary and can probably be ap- 
plied with equally good success in other forecasting problems. Spe- 
cifically, application might be very successful in attempts to predict 
final sales figures on the basis of early information, The use of 


| current data as weights in determining an average of a past per- 





formance parameter is the feature of this method, which may pro- 


vide a better current estimate of the parameter than is otherwise 
attainable. 


SAUL ROSEN 


Of the many people involved in the election forecasting effort 
special mention must be made of Mr. Thomas Jackson who pro- 
grammed and operated the UDEC computer, Professor Arvid Jacob- 
son, Director of the Wayne University Computation Laboratory, and 
Professor Robert Steadman who gave much useful advice and as- 
sistance. Professor Steadman had predicted the results with amaz- 
ing accuracy a month before the election, a feat which UDEC has not 
yet tried to duplicate. 
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Approximate Stresses in Thick-Walled Cylinders 
By Continued Fractions 
By Arvid E. Roach 


General Motors Corporation 


The purpose of this paper is to draw attention to a useful prop- 
erty of continued fractions and to show its use in an actual problem. 
An ordinary continued fraction is a fraction of the form 


1 
A = a, + i 
a, + i 
a, + 
_— | 


which for convenience may be written as 


1 1 1 
+ — _— — 

a+a,+a,+... 

Approximations to A, obtained by discarding remainders, are 
called convergents of the continued fraction. If the remainder 
—_—s -f 


>- = = be discarded, the first convergent ap is obtained; 
a + @2 + ag +... 


1 ; 
if the remainder = = be discarded, the second convergent 
2? ST ces 
+1 
ee is obtained; and so forth. 
1 


The first convergent is smaller than A, the second greater, the 
third smaller, and so on, successive values of the convergents nar- 
rowing down upon and approaching A. 

The convergents are all rational numbers, that is, they are ordi- 
nary fractions with integers dividing integers. If Pj be the numer- 
ators and Q; the denominators, then 


Po - a Pi = 4a, +1 P, = 404, 8, + aot & 
Q, °* Q, a, * Q, a,a,+1 iii 
whose numerators and denominators are, more explicitly, 
P, - ao ’ Q, = 1 , 
P, = aa, + 1, Q, 7” ee 
P, = a2)8,a,+ agt a, Q, = a@,a,+ 1 
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A systematic procedure for computing the numerators and de- 
nominators of convergents without requiring the simplification of | 
interminable complex fractions is given in the Appendix. 

It may be shown (reference 1) that any convergent é is a better 

n 


approximation to a number A than is any other rational fraction = 


whose denominator Q is less than Q,. 


This property of continued fractions is of interest to the engi- | : 
neer or designer who must carry out a large number of repetitive 
approximate computations. 

As an example, we may derive approximate equations for the 
stresses in the wall of a hydraulic cylinder. The correct stress 1 
equations are well known (reference 2): 

ail 2 2,2 : 
_ pia - b_, 1 ab (po-pi) : | 
0. b’-a r b’-a’ 
is th 
large 
2 2 24.2 
ia -pob 1 ab (po-pj) F 
o9= . —-"— ’ 
° b _ r on mame place 
exam 
in which o, is the radial stress, og is the tangential stress, a is the 
inner radius, b is the outer radius, p; is the inside pressure, and p, 
is the outside pressure. The ends of the cylinder are assumed to be 
unrestrained, is 
For the case in which p,= 0, 
o. = pia’ (1 - b’/r*)/(b* - a*) , whic! 
(1) , 
- 2 2 /y? / 2 2 (o 
Og = pjya (1+ b’/r’)/(b' - a’) , 
which shows that o, is always compressive, while og is always ten- re 
sile. The maximum tensile stress occurs at the inside surface te 
: tatior 
the cylinder wall (when r =a), or Si 
- a 2) /f,2 2 
(Fo) ax = p; (b +a )/ (b -a ). (9) 
Replacing a by R and b by R+t, one obtains (4) 
a 2 42 ’ 
(99) nax * Pil2R’ + 2Rt+ t*)/(2Rt+ t*) , (og) 





which may be converted into a continued fraction as follows: 


s the ' 


nd p, 
to be 


5 ten- 


ice of 
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2Rt +t’) 2R* + 2Rt+t? (ap = R/t +1/2 
2R’ + 2Rt + t’/2 
t’/2) 2Rt+ t?(a, =4(R/t + 1/2) 
2Rt+ t’ 


Hence, 


(09) ax 7 Pil(R/t + 1/2) + 1/4(R/t + 1/2)). 


1 
The first convergent is = +93 hence, the equation 


(g p,[R/t + 1/2] 


@) max 
is the best rational approximation formula having denominator no 
larger than 2t. 

For the stress at any other point in the wall of the cylinder, re- 
place r in Eq (1) by the appropriate value in terms of Rand t. For 
example, the tensile stress midway through the cylinder wall 


(79) = pj(8R* + 12R°t + SR*t’)/(8R*t + 12R’t? + ERE, t*) 


mid 


which in continued fraction form is 


(Ry 1 ote : 1 ] 
Pil't *-8(R/t + 172) + -i(R/t + 1/2) + 8(R/t + 1/72))" 





Although this form may not appear to be notably simpler than the 
conventional form, it is actually much more convenient for compu- 
tation, and this is its principal advantage. 

Successive approximations for the stress "e) are: 


id 
0) nia = R/t (first-convergent approximation); 
9) ‘ = [R/t-t/(8R+t)] (second-convergent approximation); 
1 


we)... * [R/t -(2Rt + t*)/(16R* + 16Rt+ 8t?)] 
(third-convergent approximation); 
: [R/t -(R*t+ Rt’)/(8R* + 12R°t + 6Rt® + t°)) 
(true theoretical stress). 
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Note that the first-convergent approximation gives results that 
are sufficiently accurate for most engineering computations and that } 
always err on the “safe” side. 

The first-convergent approximations for the stress-pressure 
ratio at the inner surface of the wall and midway through the wall 
differ by 5. The second-convergent approximations differ by 


1+ 5t/4R 
2+ t/2R’ 


which approaches ¢ as t is made smaller. 


APPENDIX 


The numerator and denominator of a convergent of a continued | 
fraction can be determined from the following formulas when the 
numerators and denominators of the two immediately preceding 
convergents are known: 


P, = a, Pa.) + Pa-2 » 
Qn ; a, Qn-! + Q,- 


Actually these formulas enable one to determine all the conver- 
gents of a continued fraction if the numerators and denominators of 
any two successive convergents are known. Two successive con- 
vergents P- ana B= 

= Q-2 


same for all continued fractions. Thus, for m = 0 





will always be known, because they are the | 


P, = a Pei + P-z = ao, 
Q, = a Qi + Q2 = 1 


For n-1, 
P, = a, P, + P.i 


Qa =3aQ@+Qi =a, 


" 
& 
3 
& 
+ 
me 


from which 


Pa * 2 and Q., = 0 
P_., = 0 and PB, = 1, 
or 
P.. 0 P 1 
2@-=e a~-.: . . . 
Q-, i and o-, 0 for all continued fractions. 
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The Life History of an Equation 


By Gerald M. Rasswejler and Carl E. Bleil 
General Motors Corporation 


This is a non-mathematical discussion of one mathematical 
equation (reference 1) 


I = KF*N’m (L. P.)e"™ A (6) V 


We have chosen this equation for a number of reasons. First, it 
is reasonably esoteric. Second, it has immediate importance in our 
laboratory at the present time because it is the foundation upon 
which we build an accurate x-ray diffraction method for measuring 
the phases in hardened steel. Third, and perhaps most important, it 
is intimately associated withthe development of a whole field of sci- 
ence, x-ray diffraction. In this development we can see the interplay 
so typical of the growth of science, the interplay between theory and 
experiment, between intuitive thought and mathematics, and between 
fundamental research and practical application. 

Whatever we call ourselves, whether it be scientist, mathemati- 
cian, or engineer, we frequently have occasion to dip into a handbook 
and pull out an equation which might look something like this one. 
We seldom give thought to its source. It is worthwhile occasionally 
to go back to examine where such equations come from, for it gives 
each of us a better understanding of what part we are playing in the 
continuous development of science. 

Before we discuss the equation itself we will introduce our ap- 
plication of it and the experimental nature of the measurements in- 
volved. This we hope will lend an air of reality to the discussion of 
the equation itself and will illustrate the type of work with which our 
Applied Science Group is concerned. We will then consider in gen- 
eral terms the broad significance of this equation and where it 
comes from — its life history. if you please. 

Figure 1 is a picture of a piece of steel taken with an electron 
microscope. The magnification here is about ten thousand times. 
This picture shows the three phases which exist in this sample, 
austenite, martensite and cementite. The angular particles are aus- 
tenite. The darker matrix is martensite. The white round particles 
are cementite, or iron carbide. The properties of the steel are ina 
large measure dependent upon the amount and distribution of these 
phases. Under some circumstances, austenite may transform at 
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normal ambient temperatures to martensite during or after fabrica- 
tion of a steel part. Since the specific volumes of austenite and 
martensite are different, this transformation may change the dimen- 
sions of the part and/or introduce internal stresses. It is apparent, 
then, that the ability to measure the relative amounts of these 
phases is basic in the metallurgy of steel. 

The classical method of determining the amount of the phases in 
a piece of steel is with the optical microscope. Under some cir- 
cumstances, the three phases can readily be seen, and the relative 
amounts measured. If, however, the individual crystal grains are 
very small, the optical microscope may fail to resolve many of the 
grains. In this case the electron microscope, which has a great ad- 
vantage in resolving power, may be brought to bear. You can see 
that the relative amounts of these materials could be obtained by 
the rather tedious procedure of measuring relative areas using 





IL THE LIFE HISTORY OF AN EQUATION 59 


micrographs of this sort. The procedure is not only time-consum- 
ing, but suffers from several rather serious limitations: (1) Only a 
very small portion of the total sample is examined here. (2) Unless 
very careful etching techniques are used, the relative areas exposed 
may not represent the relative amounts of the phases even in this 
smail area. (3) The phases may be partially transformed by the 
sample preparation. (4) The test generally destroys the part. 

Now it happens that x-ray diffraction methods are ideally suited 
to the measurement of the relative amounts of these phases. This is 
because the three phases have different crystal structures and x-ray 
diffraction is primarily a method of studying crystal structure. 
However, until recently, x-ray diffraction has not been seriously 
applied to this problem because of experimental and theoretical dif- 
ficulties. During the last few years, our laboratory has carried out 

} an intensive development of apparatus and techniques for the quanti- 

tative measurement of the phases in steel. This work has been 

outstanding and a great contribution in the application of x-rays not 
only to this probiem but to the general field of quantitative phase 
studies in all metals. Let’s see what is involved in applying x-rays 
to determining the amount of the three phases ina piece of steel. 

For the moment, let’s not worry too much about what x-ray diffrac- 

tion is. Let’s dive right into the experimental x-ray measurement. 
Figure 2 shows a schematic drawing of how the x-ray diffraction 

pattern is photographed. The x-ray beam strikes a small specimen 
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at “s” anda section of the diffraction pattern is photographed on a 
circle of film. When this film is flattened out, we obtain a series of 
lines as shown in this figure. To take pictures of this type which 
would be applicable to the determination of the phases in steel, we 
have developed and built special x-ray equipment. 


X-RAY DIFFRACTION PATTERNS OF 
AUSTENITE, MARTENSITE, AND CEMENTITE 


AUSTENITE [ 


MARTENSITE (TEMPERED) 





Figure 3 


— | 


Figure 3 shows the type of x-ray diffraction patterns which one 
obtains. Here we can compare the patterns of austenite, martensite, 


and cementite. For those who are metallurgically inclined, the fer- 
rite pattern is very similar to martensite, and in our discussion we 


are including the ferrite with the martensite. We might say that we 


have here some x-ray fingerprints of these phases. If all three 
phases exist simultaneously, the diffraction pattern will contain lines 
from all of the individual patterns. Knowing the wave length of the 
x-rays and the geometry of the apparatus, it is very easy to identify 
which of these phases are present by the positions of the lines. But 
our problem here is to measure the amount of each phase. This we 
do by taking into account the density of the diffraction lines. To ob- 
tain the densities we use a microphotometer which draws a plot of 
position versus density as shown in Figure 4. Each of these peaks 
represents a line on the diffraction pattern, the austenite marked 
“A.” martensite marked *“M,” and cementite marked “C.” So now 
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we can make a quantitative measurement of the density of these 
lines on the photograph and from this we can obtain the relative in- 
tensities of the diffracted x-ray beams. We know, of course, that 
the more austenite present, the greater the intensity and the higher 
the austenite peak on this curve. But exactly what is the quantitative 
relationship between this intensity as measured from one of these 
peaks and the amount of austenite in the sample? 


TYPICAL MICROPHOTOMETER CHART 
AND 
FILM FROM WHICH IT WAS PREPARED 
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Figure 4 


This relationship is expressed by the equation I=-GV where I is 
the relative intensity of one line in the pattern and V the relative 
volume of the corresponding phase. If we know the value of G in this 
equation for this line and phase, and if we know the corresponding 
value of relative intensity, then we can calculate the relative volume 
corresponding to V. Having calculated relative volumes for each of 
the three phases we can normalize these to percent composition, 
since we know that there are only these three phases present in our 
Sample. 

Your first in.oression might be that we could establish the value 
of G by experiments on samples of known composition. That is, we 
might measure values of I for our three phases in a series of 
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samples where values of V were known, substitute these values and 
solve for G. Unfortunately, this is not possible because there is no 
other known method for determining the relative proportions of the 
three phases in steel with the same precision we can obtain in the 
x-ray method. In other words, we have no samples in which V is 
accurately known. Therefore G must be determined theoretically 
and if our theory is wrong, our measurements will be wrong. So you 
may well imagine that various members of our staff have given very 
close attention to the theoretical derivation of what we call the “G 
factor.” We can rewrite this equation 


I = KF*N’m (L. P.) e?™ A (6) V, 


where F is the structure factor which lumps together the x-ray 
scattering power of the individual atoms and the phase relationship 
of the scattered x-rays, N represents the number of crystal unit 
cells per cm® of sample, m represents the multiplicity of sets of 
equivalent atomic planes contributing to a single diffracted beam, 


e~ represents the effect on intensity of thermal vibrations of the | 
atoms, (L. P.) is the abbreviation for the Lorentz and polarization | 
factor, and A (8) corrects for the absorption of the x-ray beam by | 


the crystal. 
Before we could apply this equation we asked ourselves many 
questions about it, just as you would do with an unfamiliar equation 


which you might lift from a handbook. We wanted to know what as- | 
sumptions, implicit and explicit, were made in the derivation. How | 


accurately could these complicated factors be evaluated? What was 
buried in the constant K, and how constant would it be in our experi- 
ments? In studying the theory to answer these questions we found 
the history of the equation to be very interesting and very closely 
associated with the development of x-ray diffraction. So let’s go 
back to the first glimmer of life in the equation, in fact, to the dis- 
covery of the phenomenon of x-ray diffraction. 

In 1912 the very existence of x-rays had been known for only 
seventeen years. Then Max Von Laue, a brilliant young Prussian, 
suggested using a crystal to determine whether x-rays were waves 
or particles. He made a tremendously important discovery by 
bringing together two things, x-rays and crystals, which had not been 
associated previously. Many important scientific discoveries con- 
sist of bringing together ideas, or phenomena, or materials, which 
have not previously had any important association in man’s con- 
sciousness. Von Laue discovered that x-rays and crystals can be 
used together to reveal much about each other. 

The Laue experiment is shown in Figure 5. A beam of x-rays 








TE 


} 


PO! 


was p 
obser 
film ; 
This 
waves 
tions ; 
acteri 
ing di 
Englar 
for stu 
diffrac 
concep 
ied in 
Spacin 
Wave |. 

So 
typical 


THE LIFE HISTORY OF AN EQUATION 

















63 

d 

o ) 

e 

e SCHEMATIC 

Ss 

ly OF A 

7 LAUE PATTERN 

G } 

ay 

lp j 

nit 

of 

m, 

the 

ion 

by SINGLE 

CRYSTAL 

any POLYCHROMATIC 

Lion SOURCE 

as- 

How Figure 5 
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pri- was passed through a crystal and a spatial interference pattern was 
yund observed behind the crystal. This was photographed by placing a 
sely | film in this space and a series of spots were obtained as shown. 
; go} This was interpreted by Laue as a diffraction effect in which the 
dis- waves scattered by each atom reinforce each other in certain loca- 

tions and interfere destructively in all others. The pattern is char- 

only | acteristic of the arrangement of the atoms inthe crystal. The amaz- 
sian, | ing discovery of Von Laue was seized upon by Lawrence Bragg in 
aves 


England who immediately appreciated the possibilities of the method 
y by} for studying crystals. 
been 
con- 


He soon formulated a simple concept of the 
diffraction phenomenon as shown in Figure 6 and expressed this 
concept mathematically in the Bragg equation which most of us stud- 
which} ied in sophomore physics, 2d sin 8 
con- 
an be 


nA This equation relates the 
Spacing, d, between parallel planes of atoms in the crystal to the 
wave length of the x-rays and the measured angle 6. 

So there was taking place in the x-ray diffraction field a very 


-rays} typical sequence in the growth of science: the discovery of a new 
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effect, the development of intuitive understanding of the effect, and the s 
its reduction to mathematical form. It was then possible to make tallog 
quantitative measurements to test and apply the theory. There is no decid 
important growth of science without quantitative measurements. prob] 
Although this concept by the brilliant young Bragg was a great count 
contribution to our understanding of x-ray diffraction, it was only an was | 
infant compared to our intensity equation. In fact, the intensity does schoc 
not enter Bragg’s equation at all, so we could not possibly use it for big s 
quantitative measurements of phase concentration. most 
When Lawrence Bragg first tried to apply his equation, he was inter: 
faced with the difficulty that the atomic spacing d and the x-ray Brag; 
wave length A were both unknown. Nevertheless, by comparing the Brill 
x-ray diffraction pattern of sodium chloride, potassium chloride, Laue. 
potassium bromide and potassium iodide, he was able to work out visitc 
the atomic arrangements and spacing in these simple crystals. E 
Using this new information on the structure of these crystals, Law- confe 
rence Bragg’s father, Sir William Bragg, built an x-ray spectrome- relat 
ter with which he obtained x-ray wave lengths characteristic of va- 1922 
rious elements. For this work, father and son were awarded the or ap 
Nobe! prize in 1915. Lawrence was, at that time, only 25 years old his ty 
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and he, I believe, is the youngest Nobelist on record even to the 
present day. 

Having some known x-ray wave lengths to work with, it was now 
possible to apply this equation to determine atomic spacing of some 
of the simpler structures. Despite the interruption by World War I, 
great progress was made in the field of applying x-ray diffraction to 
crystal structure studies during the next ten years or so. But the 
Bragg equation does not contain an intensity factor. By 1925 it was 
fully realized by all workers in the field that in order to make much 
further progress, it would be necessary to solve the problem of the 
relation between crystal structure and the intensity of the diffracted 
x-rays. 

So we have an example of another common occurrence in the 
growth of science. Experiment had outstripped theory. There were 
a number of scattered papers in the literature discussing the vari- 
ous problems involved in the intensity relationship, but no one had 
worked out a clean-cut theory which could be expressed mathemati- 
cally and put to practical use. Then in 1925, there occurred a his- 
toric meeting. The story of this meeting has recently been told by 
Paul P. Ewald (reference 2), who now heads the Department of 
Physics at the Polytechnic Institute of Brooklyn. It was Ewald’s 
doctors’s thesis which provoked the discussion that led to Von 
Laue’s discovery of x-ray diffraction. Ewald described the situation 
with regard to x-ray diffraction intensity in 1925 by saying, “All 
problems seem to end up ina sigh: If only we had a reliable means 
of measuring and evaluating intensities and of deriving from them 
the structure factor!” Ewald heard that Wyckoff, the famous crys- 
tallographer of the United States, was to be in Europe that year. He 
decided to set up an international meeting to discuss this intensity 
problem. He arranged for five days of discussion at his mother’s 
country home, twenty-five miles from Munich. The little local inn 
was rented, a blackboard was borrowed from the nearest country 
school and a few boxes of cigars placed on the table in his mother’s 
big studio. The meeting was particularly significant because for 
most of those who attended, it was the first post-war meeting of an 
international character. Those who were present were: W. L. 
Bragg, Darwin, and James from England; Wyckoff from the U. S. A.; 
Brillouin from France; Fokker from Holland; Waller from Sweden; 
Laue, Mark, Ott, Herzfeld, and Ewald from Germany. Occasional 
visitors and participants were Debye and Fajans. 

Ewald tells of the exciting days and the interesting events of this 
conference. Darwin, for example, had published two papers on the 
relation between structure and intensity, one in 1914 and another in 
1922. But they were so poorly written, few people understood them 
or appreciated their value. Darwin himself got so entangled between 
his two papers at this meeting, that he promised to restate them. 
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Before the conference closed, the problems which had to be 
solved in connection with the intensity relationships were more 
clearly understood, and Lawrence Bragg and many others attacked 
these problems with new enthusiasm and purpose. They rewrote 
Darwin’s theory. They also made experimental studies of the tem- 
perature effect which has to do with the effect of thermal motions of 
atoms in the crystal. These measurements supplement the theoret- 
ical work of Debye and Waller. It was this concerted effort which 
finally overcame the deadlock and gave us an understanding of the 
relationships which enter into this intensity equation. The applica- 
tion of x-ray diffraction to the determination of crystal structure 
has continued at a high rate since this time. 

Austenite and martensite are typical of the simpler crystal 
structures which were established quite early in the chain of events 
which we have been discussing. Austenite is a face-centered crys- 
talline form of iron. Carbon atoms may be dissolved in the aus- 
tenite occupying interstitial positions. Martensite is a body-cen- 
tered tetragonal structure which looks much like a cube stretched 
on one axis. In martensite, the carbon atoms may also occupy some 
of the interstitial positions and are responsible for the lattice dis- 
tortion. These are very simple structures, but as time went on x- 
ray diffraction techniques have been applied to much more complex 
structures. One important recent accomplishment was the deter- 
mination of the structure of penicillin in a three-year study com- 
pleted in 1945 by a group of British investigators led by Dorothy 
Crowfoot (reference 3). This is shown in Figure 7. The controver- 
sial portion of the structure is shown. The whirligig patterns rep- 
resent electronic density contours which were tediously computed 
from thousands of x-ray reflections. X-ray diffractionists are now 
engaged in an international race to establish structures of proteins 
by x-ray diffraction means. The mysteries of the chromosomes are 
also being unraveled. 

As you may well imagine, the computations involved in these 
structure analysis studies are very extensive. Since x-rays are dif- 
fracted from electrons and crystalline matter is periodic in charac- 
ter, Fourier analysis becomes applicable to solve complex crystal 
structure problems. Pepinsky (reference 4) recently devised an 
x-ray analog compuier (X-RAC) with a capacity of 1680 sine and co- 
sine functions to carry out what would otherwise be an astronomical 
number of independent calculations. Even with the computer, the 
fact that the phase relationships are unknown requires appreciable 
trial and error manipulation. In many cases, however, calculations 
with this machine have now been reduced to less time than is re- 
quired to obtain the x-ray diffraction data. 

So we have seen something of the history and the background of 
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this equation. All of the factors in this equation have undergone and 
are still undergoing careful study both theoretical and experimental. 
Some of these factors, like the structure and temperature factors, 
are fundamental in the vast structure analysis work which is now 
being carried on. X-ray diffraction is by far the most powerful 
method known for the study of physics of the solid state. 

Let us talk next about the equation more specifically as we use 
it. Let’s get some idea as to just how we obtain the values to plug 
into this equation in order to measure the amounts of the three 
phases, austenite, martensite, and cementite. 

First let us point out that we have a much simpler problem than 
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the determination of crystal structure which I have just been dis- 
cussing. We know the structure of the three crystal phases and 
must deduce the intensity relationship from this structure. 

Probably the most interesting factor in this equation is the struc- 
ture factor. To show where this comes from, let’s go back to our 
concept of what is involved in x-ray diffraction. We said a while ago 
that x-rays were scattered by the electrons in each of the atoms in 
the crystal. To build up the quantitative relationship between inten- 
sity of diffracted x-rays and crystal structure, let us start with the 
scattering of one electron as shown in Figure 8. Here the incident 
x-rays of amplitude A; strike an electron of charge e and mass m 
and are scattered in all directions. We are interested in what hap- 
pens at a point in space on the surface of a piece of film. The x-ray 
amplitude from one electron at a point was worked out long ago by 
J. J. Thompson. Depending upon the direction of the electric vector 
in the incident beam with respect to the plane of the paper, the am- 
plitude of the x-rays at the chosen point is expressed by one of two 
equations. If the beam is unpolarized, the amplitude will be the av- 
erage of these two equations as shown by the third equation in the 
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figure. So we can work out, if we wish to do so, the amplitude of the 
x-rays at a given point on our film as produced by the scattering 
from one electron in terms of the incident amplitude. 

Next let us see what the general approach is to obtaining the 
x-ray scattering by one atom. Figure 9 shows an atom model in 
which the electron density is indicated by the shading. The ampli- 
tude at the same point on our film is now expressed by the amplitude 
for one electron multiplied by an integral. The electron density U is 
expressed as a function of radius and the sine terms take care of the 
difference of phase between the electron scattering from different 
parts of the atom. This expression is integrated from the center of 
the atom to infinity and the result is called the “atomic scattering 
power f.” It can be evaluated for various values of sin 6/A\ and for 
various elemental atoms with different electron configurations. As 
the result of a great deal of work by many investigators, these 
atomic scattering powers can now be obtained from tables in stand- 
ard texts (reference 5). 

Having the amplitude due to scattering by one atom, we must now 
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add the effects of several atoms. This is usually done for one unit 
cell in the crystal, a unit cell containing the smallest number of 
atoms which can be used as a sort of building block for the whole 
crystal. In Figure 10 we assume that there are n atoms in the unit 
cell. These are not all in the plane of the paper. This figure, like 
several others, is out of scale since the distance to the film is, of 
course, not comparable to the spacing between atoms. The result is 
that the diffracted x-ray beams going to the point on thefilm are es- 
sentially parallel to each other. The procedure, in general, now is 
to add the amplitude of the diffracted beams from each of the atoms 
in the unit cell. We will pick one of these atoms for the origin. The 
amplitude at the point in question will be equal to the amplitude pro- 
duced by the atom at the origin plus the amplitude produced by the 
other atoms. But these latter amplitudes must be multiplied bya 
phase factor which depends upon the location of each atom with re- 
spect to the origin. This location is expressed in terms of a special 
coordinate system used for crystals. We may express the total am- 
plitude from all of the atoms in the unit cell as a summation from 
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j=1to n of the amplitudefor each atom, taking into account whether 
itis iron or carbon or some other element, multiplied by the phase 
factor which is determined by its position in this crystal coordinate 
system, the symbols h, k, 1, u, v, w. Since this amplitude for an 
atom is derived by multiplying the amplitude for an electron by an 
atomic scattering factor f, we may write the third equation, Figure 
11. The amplitude produced by a single electron, being the same for 
all different kinds of atoms, can be taken out of the summation. This 
summation for a unit cell is what we call the structure factor and 
appears in our intensity equation as F. 

To obtain the amplitude for a cubic centimeter of the crystal, we 
must introduce a factor N which represents the number of unit cells 
per cm*. In Figure 11 we show that the amplitude of the scattered 
x-rays fora em* of crystal is equal to the amplitude for one elec- 
tron times the structure factor times the number of unit cells per 
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cm*. Since the intensity is proportional to the square of the ampli- 
tude, these three factors are all squared in the intensity equation, 
As you can see, in our intensity equation F and the N appear squared. 
A portion of the A squared term goes into the constant. The portion 
containing l+cos 2 6 was introduced to take care of the polarization 
and goes into what we call the Lorenz and polarization factor. Ac- 
tually, we never evaluate this constant K in our work because we use 
relative values of intensity and volume. 

There are some other factors in this equation: m, which we are 
not going to discuss because we would have to introduce a number of 
new concepts; A (6), which may be considered as a correction factor 
for the absorption of x-rays and which may be calculated from ab- 
sorption coefficients; and e~™”’, which is the temperature factor to 
correct for the fact that the atoms in the crystal lattice are not at 
rest. The temperature factor has been worked out theoretically by 
Debye and Waller as we have already mentioned in discussing the 
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1925 meeting. However, it is not yet well established since experi- 
ments on the temperature factor do not agree very well with theory. 
Here is a case which is very common in applied science, where we 
must make the best use of the available information without waiting 
for complete theory, keeping in mind the limitations and the conse- 
quences. 

An example of the application of this measurement technique is 
shown in Figure 12. Here we have studied the effect of a sanding 
operation on the amount of retained austenite near the surface. This 
type of application of the intensity relationship is only a minute part 
of its overall utility. The relation between intensity and crystal 
structure, which may be expressed in various forms, has wide and 
important applications. 

Wrapped up in the mathematical shorthand of this equation are 
theories and experiments accumulated from 1912 to 1955. So, when- 
ever we lift one of these equations from a handbook we tap a vast 
reservoir of scientific knowledge. This knowledge results from the 
cooperative efforts of workers with many different skills and points 
of view. It is where the many points of view meet, the middle ground 
for theoretical and experimental science, that the mathematician and 
experimentalist should make their greatest effort. It is here where 
each finds the most difficulty and consequently where maximum co- 
operation is necessary. It is just this sort of cooperation which 
forms the underlying purpose of the Industrial Mathematics Society. 
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The Effect of Surface Residual 
Stresses on Indentation Hardness 


By Leonard G. Johnson 
General Motors Corporation 


1. DEFINITION OF HARDNESS 


Hardness is a concept concerning the definition of which there is 
no general agreement. What are its dimensions in terms of the 
basic units of length, force, and time? Every so often a private in- 
vestigator comes up with a proposed definition which best fits into 
his mental scheme of things. In fact, there are so many definitions 
of hardness that toconsider them in detail would bea lengthy project 
in itself. Without discussing any previous definitions, let us con- 
sider the following one: 

The hardness of a material under an indenter is the amount of 
biaxial effective stress reserve contained in the material between a 
state of uniform hydrostatic pressure and a state of plastic yielding. 
In other words, hardness is the additional biaxial effective stress 
beyond uniform hydrostatic pressure which must be supplied by the 
indenter before the material will exhibit plastic yielding. According 
to this definition, hardness has the same dimensions as stress, i.e., 
LBS. x IN.~* 


2. MATHEMATICAL ANALYSIS OF THE DEFINITION 


Let us assume that the surface of a specimen of material, whose 
hardness is to be measured by means of an indenter, lies in the xy 
plane of a Cartesian frame of reference, the z direction being per- 
pendicular to the given surface. Furthermore, suppose the x, y, and 
z directions are the directions of the principal stresses. 


Let Ss. = residual stress in the x direction 
(longitudinal residual stress) 
Let Ss. = residual stress in the y direction 


(transverse residual stress) 


In a state of uniform hydrostatic pressure, the principal stress 
components would be 
\ Cs * =e 


(1) < Gye * © 


Gy © = Pe ( 
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Due to the residual stresses, these become (assuming an elastic 


state) / 
\ Oo, : ~ Po 


%, = =p, * &, & 
} Jy 7 Py * Sy \ 
However, even when the principal stresses are equal to the values 
in (2), the material may not have become plastic as yet. Therefore, 


by means of the indenter we must supply additional stresses in the x 
and y directions to make the material yield. 


(2) 


Let A 


x 
Let A, 


Then the principal stress components are 


Tees, 
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the additional x stress 
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(3) 


the additional y stress 


% .* 4° p>. oS. 


Ay - Po + S, 


te 
Q 

< ” 
i} 


Under the stresses(4) the material is supposed to be plastic. Hence, 
we may assume that the von Mises yield condition holds, i.e., 


2 2 a oe 2 
(o, o o,) + (o, - o,) + (o, - G,) Sa 
where Y is the yield stress of the material in simple tension. Now, 


6,- % = (a, + Ss.) - (ay + Ss.) | 
(5) ¢ O% - 0, = A, +8, 


Wy - 0, = AY+ 8, 


Therefore, the von Mises yield condition becomes 


(6) (Ax + XP = (Ag + S.)Ay + Sy) + (Ay + s,)° = ¥° 


Assume thatA, = @ A, . Then (6) becomes 

(Ax + &)* - (Ay + SMaa, + S,) + (aa, + s,)* re , 
(7) or (1-a+a?)ay + (2-a)A,S, + (2a - IA,S, 
+ 8S, - SS, +8) = Y’ 
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stic According to our definition, the hardness of the material is the ad- 
ditional effective biaxial stress necessary to produce yielding, i.e., 
the quantity 


US err errs 


y 


= = A a , 2 
ues Let H = hardness = A, \/1- a+ a’. 
re, 
le X a 


x eS 3 
J/1-a+a?* 


Hence, (7) becomes 


(8) H? + 2- a HS, 2a - J) 


N1-a+ a’ Vilie @ + @* 


+ S, - 8,S,+S),= Y¥*. 


if H, Sx, and S, are treated as three Cartesian coordinates, then (8) 
is the equation of a right elliptic cylinder whose shape and orien- 
ce tation are functions of the number a, 


The Sx Sy projection of the cross-section of (8) formed by the 
plane H=-S, V1l-a+a’ is 


(1 - a + a’) S. + (a@ - 2) Si + (1 - 2a) &S, 
+ S, - SS, +S) = Y’ 


low, 


(9) or a*S\,- 2a@S8,S, + 8S) = Y’ 


i] 
< 


or (aS, - Sy)’ 


or of, - 4, «2 


y 


Hence, the axis of the cylinder (8) is located along the line 


Af 2 
(10) Be sAf-at ee, os 2 Bt, 


x a 7 * 


3. THE SPECIAL CASE OF A SYMMETRICAL INDENTER 


In case the indenter used in determining hardness is sym- 
metrical, it follows that it supplies equal stresses in the x and y 
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directions, i.e., Ax = Ay. Therefore, for a symmetrical indenter 
a = 1, and (8) becomes 


H*+ HS, + HSy + SX - S,Sy + S} = Y’, 


which is a right circular cylinder whose axis is located along the 
line H = - Sx = - Sy; this means that the axis of the cylinder makes 
equal angles with the positive H axis, the negative S, axis, and the 
negative Sy axis. It is therefore identical to the von Mises yield 
cylinder except that the signs in the x and y directions are re- 
versed. We can, thus, predict the effect of residual stresses S, and 
Sy on the hardness as measured by a symmetrical indenter simply 
by examining a model of the von Mises yield cylinder. 


4. CONCLUSION 


Hardness, as we have defined it, is functionally related to the 
normal residual stresses S, and S, in accordance with the second 
degree equation (8). If we solve (8) for the hardness H we obtain 


(11) H = (a - 2)Sx + (1 - 2a)Sy +4 y? = ( aSx ° Sy ) 


2V1- a+ a 2V1l- a+ a’ 








For a symmetrical indenter, i.e., for @ = 1, this becomes 


_ . (Sx_+ Sy 2 Sx - Sy\’ 
(12) w= - 58) 4/e- (eS). 





Let us put Sx M Sy = Sm = mean normal residual stress, 
~ * 
_ Sx = Sy = 8. 


Then (12) becomes 
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Approximating Functions for Digital Computers 


By Lyle R. Langdon 
Wayne University 


e 
S 
e 1, INTRODUCTION. 
d 


- The approximation of functions is not new; it has been in use 
d since man first began tothink in terms of numbers. All are familiar 


y with the approximation of 7 by as used in the grade schools; 


other examples are the terms several, many, etc. Originally, ap- 
proximation was necessary because of the lack of sufficient knowl- 
edge. Today, however, approximation is used for an entirely differ- 
e ent reason: to save time and, therefore, to save money. There 
id still is a great lack of knowledge regarding certain functions, which, 
therefore, must be approximated; for example, volume as a function 
of pressure encountered in the use of the gas tables. 
- At the Wayne University Computation Laboratory we have a 
| fairly slow machine (the UDEC), and it is imperative that we have 
good approximating functions in order that we may make a fair 
showing against the much faster machines. The faster machines 
are more costly to operate and since our approximations allow us to 
cut our time and cost it is reasonable toassume that they will do the 
same for those institutions which possess the faster machines. 

Not only must the approximating function be short and economi- 
cal of machine time; it must also have an error not greater than the 
last decimal position. In addition to the normal roundoff error there 
are three truncation errors to be considered: 


1) Error due to truncating the series. 
2) Error due to truncating the coefficients. 
3) Error due to truncating the powers of the variable. 


The last two errors are considered first when forming an approxi- 
mating function. Obviously, it does not make sense to have an ap- 
proximating function with the truncation error in the eighth decimal 
position if the other truncation errors are in the sixth or seventh 
decimal position. The error due to truncating the various powers 
may be held to a minimum if the maximum absolute value of the 
variable is kept less than 1. This may be accomplished as in the 
1 1 


case of sin x by using sin > for the range 7) to 9 OF by a change of 


variable such as — for 0< x< «. The error due to truncating 
oO 
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the coefficients is rendered null in several of the approximating 
functions by using exact coefficients. When the above two conditions 
have been met then the error due to truncating the series is con- 
sidered. 

If the series is rapidly convergent the truncation error may be 
reduced by the use of Chebyshev polynomials. However, if the se- 
ries is not rapidly convergent, one of several methods may be used, 
The first step is to analyze the function, Will a change of variable 
speed up the rate of convergence (cf. tan™x)? Is the function in the 
best possible form? It may be that a different form of the function 
will converge more rapidly ¢/. sin~x). Too much emphasis can not 
be put on the suggestion to analyze the function thoroughly. 

If the function converges very slowly I have found that a rational 
expression will quite often approximate the function satisfactorily, 
Pade approximants and the modified Taylor expansion are two 
means of obtaining such expressions. A thorough discussion of 
Pade functions can be found in The Analytic Theory of Continued 
Functions, by H. 8S. Wall and Die Lehre von den Kettenbruecken, by 
O. Perron. 


A Pade function is defined as ee where g(x) is a polynomial of 


degree p and h(x) is a polynomial of degree q such that, for f(x) an 
infinite expansion in the variable x, en will exactly equal the ex- 


pansion of f(x) up to and including the term xP*4, The error result- 
ing from truncating the series at xP*4 and using a Pade approxi- 
mant will be found to be much less than the error obtained by using 
the truncated series itself. Many analysts prefer to approximate 
functions by continued fractions. I have not found this method suf- 
ficiently economical of time; however, for those who prefer this 
method, it is recommended that Padé approximants be used rather 
than the continued fractions themselves. These approximants may 
be easily obtained by means of recursion relations, and, by testing 
successive approximants, the function may be approximated as 
closely as desired. Not much seems to have been done with such 
rational expressions since Frobenius, 

A very powerful tool for approximating some functions is a 
method originally attributed to Obrechkoff but which several more 
recent writers have rediscovered. This is a modification of the 
Taylor expansion and is well discussed and proven in the American 
Mathematical Monthly, vol. 56, April 1949, pages 243-247, by P. M. 
Hummel and C. L. Seebeck, Jr. They do not, however, show how it 
is derived. Since several of the approximating functions to be de- 
scribed were derived by means of this method it would be well to 
show how it may be derived and to list a few of the coefficients for 
the use of others who may wish to use the method. 
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2. MODIFIED TAYLOR EXPANSION, 


The Taylor expansion of f(x) about the point (a) differentiated n 
times and then multiplied by ie and arranged in powers of (x-a) 


is as follows: 


f(x) = tla) +4 as a), = + S(aXlx-a)? , 


3! 
ra) (x) = ee is altel . — re i 
ee t"(x) = Sone ‘ i eatad 
aw t"(x) = mag ail oe 


etc. 


The Taylor expansion may be approximated by stopping with the nth 
term or nt® column and the mt row. By so doing there will be m-1 
equations which will allow the m-1 following terms or columns to 
be eliminated. If m=n the work becomes much simplified, for one 


need only multiply the mt row by 1, the (m-1)th row by - er), the 
(m-2) row by im pS ald the (m-3)th row by - eee, tc. 


If we add and rearrange terms we obtain the following expressions: 


u 


n=m=1; f(x) = f(a) + s{t"(a) + £'(x) |(x-a) ; 


n=m=2: f(x) 


f(a) + s(t"(a) + £’(x) |(x-a) 
+ gel fila) - f(x) (xe-aP 


etc. 
The first n terms of the modified Taylor expansion for m=n are: 


n k k+l¢ (k) k 
(1) f(x) = f(a) + D a, (f(a) + (-1) "f(x J (S82) oR 
1 k! 2 
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The a, may be obtained from the recursion relation T 
WwW. 
_ 2(n+1- _- ; . 
k= (2nel-k) 2-15 1 = I 
It will be observed that lim a, = a,_, . 
N-> co 
Therefore Tl 


pl 
oc ak 
lim f(x) = f(a) + S[ fa) + (-1)**'4"*)¢xy } (259) 
2 
n->- 1 
It should be expected, also, that the approximating function would 
consist of two parts, each part being an approximation of (S$) 
The first term of the remaining series which is neglected gives 
a measure of the error involved in approximating the Taylor series 
bv this method, and for the case m=n is: 
Wi 
2n+1,(2n+1 
(-1)%(nt)*(xeay?* ENE ay_ is 
(2n) I nei) mat 
ar 
The remainder term for the more general case m # n is: 
ar 
(-1)®mIni(x-ay™*P*tplmtnt gy 
(m+n) ! (m+n+1)! 
Setting f(x) = e* for n=4 and a=0 it is found that 
(2) 
3 x’ 1x 1x 1x” 
(l+50 +a) +lg +a ) 
oft te 28 1680 2  &e +R 
3 x’ 1 x* 1x 1x” 
(+95 * tea) -'s * aa? 
The first bracketed expression is a very poor approximation to 
cosh 3 and the second bracketed expression is a very poor approxi- (3) 
mation to sinh 4 The resulting approximation for e*, however, is 


2 


exceedingly good, the error being in the eighth significant digit. 
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The approximation of e* obtained by this method may, therefore, be 
written 


cosh 5 + sinh 
cosh = - sinh 
2 
The bar indicates that the sub-approximation is very poor. Re- 
placing x by ix and by -ix and properly combining terms: 


2sin ~ cos ~ cos ~ - sin’ ~ 
; 2 2 2 2 
Ssinx = _ =a" cos x =— . —, 
cos ~ + sin’ * cos? ~ + sin? * 
2 2 : 2 2 
2sin > cos x 
2 2 
tan x = 
cos’ 2 sin’ = 
2 2 


Where more than one of the above functions are to be calculated it 
is recommended that they be calculated by this method, since no 
normalizing is required and after calculating one function the others 
are obtained at little extra cost. 


The approximating functions obtained by means of formula (1) 
are: 


3 [ 3 
x k x k x 
2= De-1)*b,,(=) (-1)"a,,(s) 
(2) sin x Ne GM Patg? LGM ats) J R, -1<x<+0; 
2 | 2 
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3 2k| [3 2k 
|2% 2(-1)*b2ud 3) | 2-1)kaayt 2) | 
0 ~0 —+ R, -1<¢x<+; 


3 “i 3 2k|2 
- k x te x sy_1)* x 
[2-0 azn 3) E _ I) balg) | 





3 2k 3 2k 
eaatzi"| + [3 Ebeu3) 
(5) e* =z 0 0 +R, -47¢x<+7; 
3 * r 2k | 
= tet, oe 
LR ant) }3 >baug) | 
with 
a, = 0.864864, b, = 1.297296, 
a, = 0.898128, b, = 0.37422, 
a, = 0.10206, b, = 0.0183708, 
a, = 0.0020412, b, = 0.00010935 . 


The maximum absolute values of the errors in using the ap- 
proximations in (2), (3), (4) and (5) are respectively 0.00000 0006, 
0.00000 0001, 0.00000 007, and one unit in the ninth significant digit. 


3. THE FUNCTION x?, 


The method to be outlined is valid for all real values of p, 
whether rational or irrational, and is very economical of time com- 
pared to the method which involves the use of logarithmic and expo- 
nential functions. The error can be made as small as desired. For 
XP write X™*4 where m is the integer part and |q| < 0.5. The inte- 
ger part is handled in the normal manner. Substitute x49 in formula 
(1) of the modified Taylor expansion and rearrange in terms of “., 
which we denote by y. It will be found that further simplification is 
possible if, as in the case of the trigonometric functions, the re- 
sulting expression is subdivided into two separate functions which 
we denote by A and qyB. Then 


. [An + ayBn] 
(6) ” [An - qyBn| *  ’ 





wh 
sat 


(7) 
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where n is the number of terms used in formula (1). Both A, and B,, 
satisfy the following recursion relation: 


(7) an eS yc 


2 
Starting with A, =1,A,=1+ (aoy , B, = B, = 1, x9 may be ap- 


proximated to any degree of accuracy desired. When q = 0.5 it will 
be observed that the recursion relation simplifies to: 


Cn = Cn-1 - Cn-2- 


If qy is grouped as a new variable the resulting expression may be 
written in terms of binomial coefficients, as for example: 


0-5 _ 0-5 [1+y] S$ . 0.5 {( 1-y*)+y] 
x a fi-y| x” a (1-y)-y] ? etc. 


where - 1(x-a) 


2 (x+a) 


With a = 1 the following approximation will yield the square root of x 
with a maximum error in the fifth significant digit for n=7 and 


01¢ x<¢ 10.0. 


“ = A + yB = 1 (x- 1) 
: A - yB ’ 2 (x+1) ’ 
3 3 
A = Z(-1)*a,y** , BB = 2{-1)*b2xy2k 
0 
a = 1, b, = I, 
a, = 6, b, ~ 5, 
a, = 10, b, = 6, 
a, = 4, b, = 1 


If the result of this approximation is used once in the Newton- 
Raphson equation the resulting error will be outside the range of the 
average machine. 
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; ; : X-a a ‘ 
Some time is lost in forming —— but this is necessary in order 


to avoid exceed capacity and loss of digits out of the right side. Any 
ene of the approximants may be used as an iteration method for q 
equalling a reciprocal of an integer (a # 1); the x9 becomes the 
new a%. Obviously, this is not the case for q, a number other than a 
reciprocal of an integer; then a must be equal to 1. The recursion 
formula will be used up to and including the value of n for which the 
error is permissible. 


4. APPLICATION OF CHEBYSHEV POLYNOMIALS. 


I shall make no attempt to enter into the theory of Chebyshev 
polynomials; the reader will find an excellent treatment of them in 
“Trigonometric Interpolation of Empirical and Analytical Func- 
tions,” by C. Lanczos and in “Tables of Chebyshev Polynomials,” by 
the National Bureau of Standards in the Journal of Mathematics and 
Physics, vol. 17, 1938, Some explanation, however, is necessary in 
order that the reader may further truncate some of the approxi- 
mating functions to be given, where the conditions of his problem 
will allow a greater amount of error. 

The Chebyshev polynomial of degree k may be written: 

Tx) 
gk- 1 


= x™ + a, _oxk-2 + ay _gxk-4 +... + a 


, 
or we may write 


Tx(x) 


k-4 
tus * & 9 9k- | 


(8) x* = -(a,_px*-2 + ay _4x 


The a; may be obtained from the tables mentioned. The function 
f(x) = bo + byx + box” +.. + b,x" may be truncated or reduced to a 
polynomial of a lower degree by replacing x™ by equation (8) if it is 
found that the error committed by neglecting the term PaTo(x) is 
permissible. The absolute value of T,(x) will never exceed 1. This 
process may be continued until the error involved will no longer 
permit such truncation. 

An illustration of the use of the Chebyshev polynomials will be 
helpful for those who are unacquainted with them. 

Tn(x) equals cos nQ; therefore, before any attempt is made to 
use the Chebyshev polynomials, the variable must first be assigned 
the definite range Tl or a transformation must be made so that the 
range will be 11, since cos n@ ranges from -1 to +1. For our 





APF 


exal 
to + 
mat 
toe 


Us: 














APPROXIMATING FUNCTIONS FOR DIGITAL COMPUTERS 87 
example we shall use sin x and stipulate that x will range from -1 
to +1. An approximation polynomial is desired which will approxi- 
mate sin x between -1 and +1 with a maximum truncation error not 


to exceed 0.00000 005. The expansion is 


x? x? ; » r x? : x! + 
5040 362880 399 16800 °** 





The 11th degree term does not exceed 0.00000 005, so we may ne- 

x! 
399 16800 
exercised at this point; we must be sure that the summation of all 
terms following does not exceed the permissible error. For the 
function sin x we know that the maximum error will be less than the 
first term neglected. From the Tables of Chebyshev polynomials 
we have 


glect since it actually is 0.00000 00250- - . Care must be 


o _ 576x" 432x° 120% 9x T,(x) 
x _—_—— 


" 356 ~ 256 * 356 ~ 256° 356 ° 
x? = Lax’ _ 56x | Ix, TAX) 
64 64 64 64 


Using the equation for x” we have 


x” . x” x’ 
initia i it 
‘ail (ay 432x',120x 9x , Tdx) ) 
3 62880 \256 256 256 256 256 


After multiplying and combining terms 


928 97271x 154 82760x° 





sin X = 99897280 ~ 928 97280 
7 73712x° 17856x" T(x) 
* 92897280 “92897280 *928 97280 ° 
The maximum error committed by replacing x” is FAX) or 
Se 928 97280 


0.00000 00108; therefore, we may neglect this term and write sin x 

= 0.99997 88799x-0.16649 71892x*+0.00799 23115 + R, where |R| < 

0.00000 0011. Using the same procedure we might continue by using 
1, ; 17856 TAx) _ 

the equation for x’; but, we find that 928 9728064. 0.00000 30034 


for T,(x) equalling 1. This is an error greater than we can tolerate; 
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therefore, we go no further. The approximating polynomials which 
follow were formed in this manner. They have a very small trun- 
cation error and may be used as they are or they be truncated 
further as indicated above. 


Sin 5 = Zeon ee +R, -1<x<+l1, Cos = = Ze a_n + R, 
c, = 1.57079 62899, co = 0.99999 99998, 
c, = -0.64596 33587, c, = -1.23370 05336, 
cs = 0.07968 84770, c, = 0.25366 93147, 
(8-1). = 9.00467 22243, '%7) © = 0.02086 26564, 
co = 0.00015 08196, cs = 0.00091 76703, 
IR! < 0.00000 0004 . Cio = -0.00002 37888, 


IR! < 0.00000 0007 


The function logex can be evaluated in a similar manner. The 
method to be outlined will permit the evaluation of the logarithm 
with a minimum amount of coding and computing time: 


6 
(8.3) log.x = loga + Ecz,,,;y°"*!' +R, 
o 


c, = 2.00000 00366, 
c, = 0.66666 17100, 
c, = 0.40019 30326, 
c, = 0.28243 35712, 
Cy = 0.25034 10930, 
C, = 0.05722 83265, 
C3 = 0.41059 70438, 
IR! < 0.00000 00003 


In this series y = — where a is a function of the square root of 10. 


Ht O10<x< 10, a = 2. 
Hf 100<x< 100,a= V10. 
If 10.00 < x < 100.0 , a = 10/10 , etc. 


(10) te 


AS xn 
of 0.41 

Sin 
expans 
aforen 








10. 
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5. THE FUNCTION ARCTANGENT X, 


The method to be outlined will provide a very economical means 
of computing the values of the arctangent; moreover, the truncation 
error is exceedingly small. We shall start with the well known re- 
-1 X +X, 


lation tan~'x + tan™'x, = tan . 
1 - xx, 


An effort will be made to find 


x+xX,, a : 
the value of x, for which 1 — is a minimum in absolute value and 
= 1 


will have the same absolute value at x=0 as at x=1. Once this is 
accomplished the series expansion from the arctangent may be 
truncated by means of the Chebyshev polynomials. 


We shall let y represent the expression — = 

™ 1 

y will equal x, and for x=1 we shall arbitrarily set y equal to -x. 
1 +x, 


The result is 2 = -~-x,. The solution of the quadratic yields x, 
- 1 


; then, for x=0 


= 1+N2. Since we wish the value to be as small as possible we 


choose x, = 1 - V2 = - 0.41421 35623-- . With this choice the equation 
becomes 


_ X - 0.41421 35623 


-1 ap -1 
meee gt ey. Te aeel eees SS XS 


8 1. 
As x ranges from 0 to 1 the new variable y will range from a mini- 
mum of -0.41421 35623 to 0 and then increase to a maximum of 
0.41421 35623. 


For values of x greater than 1 the equation will still be valid if x 


1 ' a 1 , 
is changed to ~ and the equation tan™ x = 3 tan~* = is subtracted 


2 is eliminated. The result- 


from this new equation (9) so that tan™ 


ing equation is 


1 - 0.41421 35623x 


an “= ©. os 
(10) tan“'x = —j - tan”y , ¥ = [0 43421 35623 1S *S *- 


8 





As x now ranges from 1 to the variable y ranges from a maximum 
of 0.41421 35623 to 0 and then to a minimum of -0.41421 35623. 
Since the maximum absolute value of y is very small the series 
expansion for arctangent x may be economized by means of the 
aforementioned Chebyshev polynomials. The resulting equation is: 


4 
-1 


- ae 2n+1 
tan y = ZCans1¥ 
oO 


RR, 
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c, = 0.99999 99031, 
c, = -0,.33332 18453, 
c, = 0.19961 55679, 
c, = -0.13751 64984, 
c, = 0.07726 42020, 
IRI < 0.00000 0004 . 


y is the result of substituting into the equations the value of x: 


x - 0.41421 35623.) 

1 + 0.41421 35623x ’ °S*<! » 
and 

1 = 0.41421 35623x ye ew 


x + 0.41421 35623 ’ 


Since there are only five terms in the polynomial and only two 
ranges requiring but one test the degree of economization attained 
is exceedingly high. 


6. THE ARCSIN FUNCTION. 


, — , d 
The basic equation for the arcsin is: sin™'x = / 1 — If we let 


x = cos 9, then dx = - sin 9 dQ and sin™'x = - / dQ. Since ‘ 


V1 - cos*n@ tht. 


this may be inserted under the integral sign and as a result sin™ 

sin n@ de , ; 1 

S « x j d -1 Se 

S Ji - eoman Integration of the right side yields sin™x : 

sin~cos n@ + C, But cos @ =x and cos n@ is the Chebyshev poly- 
nomial of degree n; therefore, the expression may be written 


i 


x 


_ (n-1)7 1 
n 


te | 
on sin ‘ T,,(x) 


(11) sin™'x 


The equations which will be used are: 


(a) sin™x = sin~' T,(x) 0< x< 0.5 
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(b) sin™x = ; + ; sin T,(x) O5<xX< 0.866 
(c) sin™“x = = + ; sinT,(x) 0.866< x< 0.966 
(d) sin “x = ; "+ ; sin~T,(x) 0.966 <x< 1.0 


When dealing with the arcsin exercise care since we are dealing 
with a multiple valued function and must not go from one branch to 
another. No trouble will be experienced in the method being outlined 
provided that one does not go outside the ranges indicated. The 
functions T; are as follows: 


T, = x, T, = 2x*- 1, T, = 4x° - 3x, T, = 8x*~- 8x*+ 1 


For the ranges indicated, T,, T, and T, will range from -3 to +3: 
Using the series expansion for sin™ T,(x) 


wo 7 Ix’ 1°3x° = 1-3-5x" 

ed ee * ee 
where x will never be greater than 0.5 in absolute magnitude. If we 
replace x by 5 then the new variable will range from -1 to +1 and by 


let the use of Chebyshev polynomials the series may be converted to an 
approximating polynomial. When we have done this, we then trans- 


form x back to its original value ranging from -3 to +3 and 


5 


(12) sin”x = Xe,,,,x**t! +R, ~7 S25 o4 . 
b o- =. Oe 
- 
1 
“1 where c, = 0.99999 99711, 
oly- c, = 0.16666 98337, 


c, = 0.07490 14744, 
c, = 0.04593 87798, 
c, = 0.02231 69693, 
c,, = 0.04485 69846, 
< 0.00000 0003 
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placed by the value of T,(x). The expression as given is for T,(x) 
since T,(x) = x. 

The method outlined provides a very economical means of com- 
puting the arcsin up to 0.966. It is not economical to use higher de- 
grees of T,{x) in order to approach closer to 1. For this range two 
methods are available. 


Method I. 
This method uses the equation 
* 


1 - x 
2 { 2 


sin"x = 


Since (1 - x) is near zero, economization of the series by means of 
the Chebyshev polynomials is indicated. The equation then becomes 


4 
- V1i-xZoc,(1 - xk +R, 0.966<x<1.0 , 
Oo 


(13) sin™x = 


wlsa 


Co = 1.41421 35625, 
c, = 0.11785 10948, 
C2 = 0.02651 86007, 
c; = 0.00784 85583, 
C4 = 0.00304 49584, 


Due to the fact that (1 - x) is near zero the computation of the 
root will take longer and the error will be greater than by the 
method following. 


Method II. 


This method uses the negative portion of the curve of T,(x), 


sin~’ [3x - 4x*] = sin 'x = 3 sin“s . 
x is the root of the equation 4x°- 3x +x = 0. It may be obtained by 
the standard Newton-Raphson method for which a very good first! 
approximation is available since 0.4227< x< 0.5 for 0.966< 1 
< &e. 
One need not be perturbed by the zero of the denominator for 


It is understood that the x’s in the above equations are to be re- | 


\R\ < 0.00000 0003 . 





is by 
four | 
fore, 


(14) 


(15) F 
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aan x = 0.5 since the numerator tends to zero faster than the denomi- 
.) nator. Example: 


\- required sin™'(0.99999 999) . x = 0.49995 904 


, 


10 therefore, sin~'(0.99999 999) = 3 sin™*(0.49995 904) . 


7, THE ERROR FUNCTION. 


The most economical means of computing the function 


H(x) = - L* e-*'dx 
is by separating the function into four ranges. The result is a set of 
of four polynomials which involves no computation of ex” and, there- 
- fore, represents a considerable saving of computation time. 
2 Xx 2 6 2 
(14) H(x) * i+ Tres F en eR, O52 SI, 
oO 
a, = 1.12837 91582, 
a, = -0.37612 60541, 
as = 0.11283 42553, 
a, = -0.02684 84339, 
ag = 0.00517 94871, 
ai. = -0.00079 40256, 
the ais = 0.00007 64065, 
the IR! < 0.00000 00006 . 
' (15) H(x) = FLX e-*"Ax = ao + kag (2-8) +R,1<x<2, 
2 
a, = 0.96610 51464, 
a, = 0,99999 99787, 
od by a, = -0.74999 99270, 
bre a, = 0.29166 70928, 
a, = -0.04687 56199, 


-0.00781 48944, 
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a, = 0.00508 00320, 


a, = -0.00061 76364, an 
a, = -0.00015 77822, se 
a, = 0.00005 05946, shift 
k = 0.05946 51446, 
IRI < 0.00000 00002 ¢ 
T 
(16) H(x) = J Zen dx = a,+ kea,(2x-5)" BR, 34 2S87 >< 
5 1 and 
a, = 0.99959 30388, bon 
a, = 1.00001 55112, (18) 
a, = -1.24973 14425, 
a, = 0.95812 62601, (19) 
a, = -0.49611 57754, 
a, = 0.17626 77866, (20) 
a. = -0.03808 98415, (21) 
a, = 0.00289 10547, 
k = 0.00108 91421 15, rm 
IRi < 0.00000 0008 . ee 


2 ,x 6 
(17) H(x) = de 7 eX dv = ao + = ean +R, 34 x<4, 
2 
a, = 0.99999 92580, 
a, = 1.02077 67980, 
a, = -1.76324 47835, 
a, = 1.79422 19871, 
a, = -1.50194 84028, | 
a, = 1.27034 44487, 
a, = -0.55367 23595, ' 
k = 0.00000 26997 13388, 


0.00000 0007 . 
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For all practical purposes the value of the function may be con- 
sidered equal to 1 for all values of x greater than 4; the maximum 
absolute value of the error is less than 0.00000 0016. In order to 
obtain the accuracy indicated it is necessary that the constant k be 
shifted as far to the right as possible before multiplying. 


8, e*sin x AND RELATED FUNCTIONS. 
The functions e*sinx,e*cosx, /e*sinxdx and “ e*cosxdx all satisfy 
\ O° ‘ 
the differential equation y(*)+ 4y = 0. If y is set equal to = a., .x°*™ 
) 


and substituted into the differential equation four solutions will 
result. 


s0(. 4) y 4n 
(18) y, = x Tah = coshxcosx , 
19 2 Se. coshas sinhxc 
(19) y, = > (4n+1)! = coshxsinx + sinhxcosx , 
>0(. 4)" x Ante 
= >» — = 
(20) V3 > [ane2)! sinhxsinx , 
s0f_ 4)" y 4n+3 
(21) y, = 43 . = coshxsinx - sinhxcosx . 


o \4nt 


Now let each of the y; be approximated over the range - L to on 


2 2 
by a polynomial using Chebyshev polynomials for the purpose: 


y, = - sg oR , y, = eee 
oO ) 

a, = 0.99999 99999, a, = 3.14159 26536, 

a, = -1.01467 80274, a, = -0.63754 10092, 

a, = 0.01470 81434, a, = 0.00513 41140, 

a. = -0.00003 01160, a,,; = 0.00000 72798, 

IRI < 0.00000 00001 . IRI < 0.00000 00002 

Ys = Sadnsaxtt2 +R, ¥ * Saaneaxt™t3 +R, 
o o 

a, = 2.46740 10108, a, = 2.58385 63712, 


-0.16690 71308, -0.07490 79164, 
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ao = 0.00080 50348, a, = 0.00023 00264, In t 
IR| < 0.00000 001 . IRI < 0.00000 0003 


The y; are the elements of the following eight functions: 


™ 


coshs cos 5 x =y , eZ sin 5x = Hyatys) +y, , oe 
from 
PO, — 7 * 7 1 T 
cosh xsin5 x + sinh 5 Xcos 5 x = hh, e* COS5x = V2 Va) oy « a 
tage 
sinh5 xsin5x = ys , Sexsinxdx = Hv +h -¥1) both 
' posit 
bl i bd o j Ld Ld = x = 1 time 
cosh 5 xsin 5x sinh 5 Xcos 5 x ¥,,/e*cosxdx BV. t+¥ ats) daar 
tion; 
The advantages of the above method are: roots 
1) Each of the original functions has a set of terms missing; nothing 

is done to disturb this condition. 10, 1 

; 
2) When more than one of the above functions are required the Al 
others are obtained at very little extra expense. | signif 
theret 
_* 1 : by th 
3) The relation e* sin (5 + x) = e2(e*cosx), etc. allows one to ae 
compute the functions for any desired x by scaling the constant that a: 
term. Table 
run tc 
| Tate t 
9. AN INTERPOLATION EQUATION. } in the 
equati 


A computation laboratory has frequent need for an interpolation| ,, ma 
equation and the one which is about to be described is superior t0 | 4, 
the standard Lagrange equation. The number of additions and multi- Sm 
plications necessary with the latter make it very uneconomical ad} accom 
time; a great amount of storage is also required. | 

; nearly 
= (x-X, (X-X, )foH{ X-Xo)(X-X,) aif +H x-X.xX-x, Jaf, varies 
(x-x, (x-x,) Hx-x)(x-x,)a, Hx-xQ(x-x,)a, | — 


The equation is exactly satisfied at the points x,,x , and x,; however, T=400‘ 
a, and a, must be evaluated. Since this is a fourth degree equation} 
two more values f, and {, at the points x,and x, must be used. It is 1,952 a 


suggested that the points be chosen thus: | partial: 


the 
Xo Xs X XX XX ran 


(22) Let f(x) 





YON 


lation 
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multi- 
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In the case of equal intervals the coefficients of the a; are much 
simplified. Rearranged, equation (22) becomes 


(x-x,)(f-f)) 
4. * (x-x\f-f,) 2 


(x-x,)(f-f) 
~ (x-x,)(f-f,) 


Substitution of the values of x,,x,,f, and f, will give two equations 
from which to solve for a, and a». 

This rational expression interpolation equation has a much 
smaller error term than the Lagrange equation; but its big advan- 
tage is the great saving of time. In one problem the UDEC formed 
both a fourth and a sixth degree equation and tested for the error 
position and used the equation which best suited the problem. The 
time required was less than that previously necessary for a fourth 
degree Lagrange equation. This is strictly an interpolation equa- 
tion; extrapolation must not be attempted! The denominator has 
roots just outside the range. 


10, THE GAS TABLE NO.1, by KEENAN AND KAYE, 


About half of the tabulated data in the gas tables is given to five 
significant figures, the remainder being given to only four figures; 
therefore, it is not possible to exactly reproduce five digit numbers 
by the use of four digit numbers even after roundoff as the tables 
now exist. An effort must first be made to smooth the values such 
that after roundoff the tables will be reproduced. This was done for 
Table I, but owing to the lack of machine time the smoothing was not 
run to completion. The equations formulated are, therefore, accu- 
rate to only four figures. Owing to the importance of the gas tables 
in the field of heat engines it is thought advisable to include the 
equations at this time even though they are not as good as they may 
be made. Some details of the procedure will be given for the benefit 
of those who may wish to continue the work. 

Smoothing of data and the formulation of equations may best be 
accomplished if the function is first put into a form which is as 
nearly linear as possible. If f is tabulated it will be found that it 
varies only slightly from linearity, being 0.238825 at T=400° to 
0.2605 at T=2800°. The ratio me (scaled) varies from 1.89765 at 
T=400° to 1.1421 at T=2800°. The ratio (scaled) vT* varies from 
L952 at T=400° to 3.2434 at T=2800°. The functions + and ©. were 
partially smoothed and the approximating functions formulated for 
the range T=400° to 2800”. 
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(24) 


us 


nw 
OM e/OM. 


0.24733 86410, 
-2.78668 77031, 
-0.13588 75904, 
-5.03625 92978, 
1.66128 75237, 
0.50000 00000, 
-5.71356 80651, 
0.53688 35871, 
-9.03206 87779, 


4.41184 55725 . 


4h 


OMY OM. 


0.50360 91654, 
-3.60204 13117, 
- 1.60756 25738, 
-4.46070 90814, 
2.09553 00034, 
0.50000 00000, 
-3.49851 89290, 
-2.07214 43879, 
-5.28452 66199, 


1.37873 08856 . 
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- 22” 

,. = oe 
5 oil 

T = 4000 ’ 


400° < T< 2800° , 


ho = 414.507944 


’ 


95.53 < h< 732.33 





APP 
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(26) 


(27) 
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95.5301688 < h < 395.7418976 , 


4 
=: n |a 
$ 10-"* ant . : he ' 

. any” ho = 182.080895 , 
= 0.026666 99359, a) = 0.003314 19725, 
= 0.100133 5247, a, = 0.013280 65434, 
= 0.194058 417, a, = 0.027305 2641, 
= 1.021617 398, a, = 0.133409 354, 
= 0.058822 6247, a, = 0.039548 6036, 


The quotient of the two polynomials lies between 7.35 and 8.74. 


? (26) 


(27) 





3 ‘ 
hE bpy™ . _h-h, 
¥\-8 10-** h+h,’ 
= any" ho = 560.59411 54 , 
395.7418976 < h< 513.4 
= 0.17848 88660, a) = 0.02500 00000, 
= -1.27100 55585, a, = -0.17450 20952, 
= -1.44345 68435, a, = -0.23142 07406, 
= -1,79322 13828, a, = -0.24942 69209, 
3 n 
—~ 4) 2 thy _ P- Po 
= p98 10° a: 
= ayy" Po = 12,29821 93 , 


= 
~ 


0.48577 79< p < 71.72885 18 , 
= 1.28686 75130, a = 1.00000 00000, 
0.56780 65679, 0.37230 55663, 
-0.75104 40286, = -0.60632 45096, 
-0.18877 74265, -0.12366 56179, . 


rT 


Pr KP 
i 
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» ©- Be 


% " p+ Do” 


2 
= 
) 
2 
z 
9) 


71.72885 18< p< 702.000897 , 
1.40064 75900, 1.00000 00000, 


P, = 256.612982 


, 


1.08112 85249, 0.73384 35822, 
-0.05531 60092, -0.05471 70208 . 
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